Numerical solution of the Holstein polaron 
problem 



H. Fehske 1 and S. A. Trugman 2 

1 Institut fur Physik, Ernst-Moritz-Arndt-Universitat Greifswald, D-17487 
Greifswald, Germany holger.fehske@physik.uni-greifswald.de 

2 Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 
87545, U.S.A. sat@lanl.gov 

1 Introduction 

Noninteracting itinerant electrons in a solid occupy Bloch one-electron states. 
Phonons are collective vibrational excitations of the crystal lattice. The basic 
electron-phonon (EP) interaction process is the absorption or emission of a 
phonon by the electron with a simultaneous change of the electron state. From 
this it is clear that the motion of even a single electron in a deformable lattice 
constitutes a complex many-body problem, in that phonons are excited at 
various positions, with highly non-trivial dynamical correlations. 

The mutual interaction between the charge carrier and the lattice defor- 
mations may lead to the formation of a new quasiparticle, an electron dressed 
by a phonon cloud. This composite entity is called a polaron [1, 2]. Since the 
induced distortion (polarisation) of the lattice will follow the electron when it 
is moving through the crystal, one of the most important ground-state prop- 
erties of the polaron is an increased inertial mass. A polaronic quasiparticle 
is referred to as "large polaron" if the spatial extent of the phonon cloud is 
large compared to the lattice parameter. By contrast, if the lattice deforma- 
tion is basically confined to a single site, the polaron is designated as "small" . 
Of course, depending on the strength, range and retardation of the electron- 
phonon interaction, the spectral properties of a polaron will also notably differ 
from those of a normal band carrier. Since there is only one electron in the 
problem, these findings are independent of the statistics of the particle, i.e. 
we can think of any fermion or boson, such as an electron-, hole-, exciton- or 
Jahn- Teller polarons (for details see Refs. [3, 4, 5]). 

The microscopic structure of polarons is very diverse. The possible situa- 
tions are determined by the type of particle-phonon coupling [4, 6]. Systems 
characterised by optical phonons with polar long-range interactions are usually 
described by the Frohlich Hamiltonian [7, 8, 9]. If the optical phonons have 
nonpolar short-range EP interactions, Holstein's (molecular crystal) model 
applies [10, 11]. For a large class of Frohlich- and Holstein-type models it has 
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been proven that the ground-state energy of a polaron is an analytic function 
of the EP coupling parameter for all interaction strengths [12, 13, 14, 15]. The 
dimensionality of space here has no qualitative influence. In this sense a (for- 
mal) abrupt (nonanalytical) polaron transition does not exist: The standard 
phase transition concept fails to describe polaron formation. It is, instead, a 
(possibly rapid) crossover. (We mention parenthetically that in contrast to 
the ground state, the polaron first excited state may be nonanalytic in the 
EP coupling.) 

The fundamental theoretical question in the context of polaron physics 
concerns the possibility of a local lattice instability that traps the charge 
carrier upon increasing the EP coupling [1]. Such trapping is energetically 
favoured over wide-band Bloch states if the binding energy of the particle ex- 
ceeds the strain energy required to produce the trap. Since the potential itself 
depends on the carrier's state, this highly non-linear feedback phenomenon is 
called "self-trapping" [3, 4, 16, 17]. Self-trapping does not imply a breaking of 
translational invariance. In a crystal the polaron ground state is still extended 
allowing, in principle, for coherent transport although with an extremely nar- 
row band. One way to think of this is that a hypothetical self-trapped state 
can coherently tunnel with its phonon cloud to neighbouring locations, thus 
delocalising. The problem of self-trapping, i.e. the crossover from rather mo- 
bile large polarons to quasi- immobile small polarons, basically could not be 
addressed within the continuum approach. Self-trapping requires a physics 
which is related to particle and phonon dynamics on the scale of the unit 
cell [18]. On the experimental side, an increasing number of advanced mate- 
rials show polaronic effects on such short length and time scales. Self-trapped 
polarons can be found, e.g., in (non-stoichiometric) uranium dioxide, alkaline 
earth halides, II- IV- and group- IV semiconductors, organic molecular crystals, 
high-T c cuprates, charge-ordered nickelates and colossal magneto-resistance 
manganites [19, 20, 6, 21, 22, 23, 24, 25, 26]. 

As stated above, the generic model to capture such a situation is the 
Holstcin Hamiltonian, which is most simply written in real space [10]. Here 
the orbital states are identical on each site and the particle can move from site 
to site exactly as in a tight-binding model. The phonons are coupled to the 
particle at whichever site it is on. The dynamics of the lattice is treated purely 
locally with Einstein oscillators describing the intra-molecular oscillations. 

Theoretical research on the Holstein model spans over five decades. As 
yet none of the various analytical treatments, based on variational ap- 
proaches [27, 28] or on weak-coupling [29] and strong-coupling adiabatic [10] 
and non-adiabatic [30, 31] perturbation expansions, are suitable for the in- 
vestigation of the physically most interesting crossover regime where the self- 
trapping crossover of the charge carrier takes place. That is because precisely 
in this situation the characteristic electronic and phononic energy scales are 
not well separated and non-adiabatic effects become increasingly important, 
implying a breakdown of the standard Migdal approximation [32] . The Hol- 
stein polaron can be solved in infinite dimensions (D = oo) using dynamical 



Numerical solution of the Holstein polaron problem 



3 



mean-field theory [33, 34]. While this method treats the local dynamics ex- 
actly, it cannot account for the spatial correlations being of vital importance 
in finite-dimensional systems. 

In principle, quasi approximation-free numerical methods like exact diag- 
onalisation (ED) [35, 36, 37, 38, 39], quantum Monte Carlo (QMC) [40, 41, 
42, 43, 44, 45] and diagrammatic Monte Carlo [46] simulations, the global- 
local (GL) method [47] or the recently developed density-matrix renormali- 
sation group (DMRG) technique [48, 49] can overcome all these difficulties. 
Although most of these methods give reliable results in a wide range of pa- 
rameters, thereby closing the gap between the weak and strong EP coupling, 
low- and high-frequency limits, each suffers from different shortcomings. ED is 
probably the best controlled numerical method for the calculation of ground- 
and excited state properties. In practice, however, memory limitations have 
restricted brute force ED to small lattices (typically up to 20 sites). So results 
are limited to discrete momentum points. QMC can treat large system sizes 
(over 1000 sites) and provide accurate results for the thermodynamic proper- 
ties. On the other hand, the calculation of spectral properties is less reliable 
mainly because of the ill-posed analytic continuation from imaginary time. 
The GL method is basically limited to the analysis of ground-state proper- 
ties. DMRG and the recently developed dynamical DMRG [50] have proved to 
be extremely accurate for the investigation of ID EP systems, where they can 
deal with sufficiently large system sizes (e.g., 128 sites and 40 phonons). The 
determination of spectral functions (in particular of the high-energy incoher- 
ent features), however, is computationally expensive and so far there exists 
no really efficient DMRG algorithm to tackle non-trivial problems in D > 1. 

In this contribution we provide an exact numerical solution of the Holstein 
polaron problem by elaborate ED techniques, in the whole range of parame- 
ters and, at least concerning the properties of the ground state and low-lying 
excited states, in the thermodynamic limit. Combining Lanczos ED [51] with 
kernel polynomial [52, 53] and cluster perturbation [54, 55] expansion methods 
also allows the polaron's spectral and dynamical properties to be computed 
exactly. A numerical calculation is said to be exact if no approximations are 
involved aside from the restriction imposed by finite computational resources, 
the accuracy can be systematically improved with increasing computational 
effort, and actual numerical errors arc quantifiable or completely negligible. 
In most numerical approaches to many-body problems, the numerical error 
decreases as l/log(effort), where effort means either execution time or storage 
required. Thus even a large increase in computational power will not greatly 
improve the accuracy. Despite some progress by virtue of DMRG-based basis 
optimisation [56] or coherent-state variational approaches [57, 58], ED of EP 
systems remained inefficient. Recently ED methods have been developed that 
converge far more rapidly, with error ~ 1/ (effort) 9 , where 9 is a favourable 
power (0 w 3 at intermediate coupling) [59]. Thus doubling the size of the 
Hilbcrt space results in almost an extra significant figure in the energy. The 
algorithm [60, 61, 62] we will apply in the following is based on the construe- 
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tion of a variational Hilbert space on a infinite lattice and can be expanded 
in a systematic way to easily achieve greater than 1 0-digit accuracy for static 
correlation functions and 20 digits for energies, with modest computational 
resources. The increased power makes it possible to solve the Holstein polaron 
problem at continuous wave vectors in dimensions D=l, 2, 3, 4, ... . 

The paper is organised as follows: In the remaining introductory part, 
Sect. 2 presents the Holstein model and outlines the numerical methods we 
will employ for its solution. The second, main part of this paper reviews our 
numerical results for the ground-state and spectral properties of the Holstein 
polaron. The polaron's effective mass and band structure, as well as static 
electron-lattice correlations, will be analysed in Sect. 3. Section 4 is devoted 
to the investigation of the excited states of the Holstein model. The dynamics 
of polaron formation is studied in Sect. 5. Characteristic results for electron 
and phonon spectral functions will be presented in Sec. 6. The optical response 
is examined in Sec. 7. Here also finite-temperature properties such as activated 
transport will be discussed. In the third part of this paper finite-density and 
correlation effects will be addressed. First we investigate the possibility of 
bipolaron formation and discuss the many-polaron problem (Sect. 8). Second 
we comment on the interplay of strong electronic correlations and EP inter- 
action in advanced materials (Sect. 9). Some open problems are listed in the 
concluding Sect. 10. 

2 Model and methods 
2.1 Holstein Hamiltonian 

With our focus on polaron formation in systems with short-range non-polar 
EP interaction only, we consider the Holstein molecular crystal model on a 
D-dimcnsional hyper-cubic lattice, 



where c] (cj and b\ (6J are, respectively, creation (annihilation) operators 
for electrons and dispersionless optical phonons on site i, and n i = c-Cj is 
the corresponding particle number operator. The parameters of the model are 
the nearest-neighbour hopping integral t, the EP coupling strength g, and the 
phonon frequency ojq. The parameters t, Wo, and g all have units of energy, 
and can be used to form two independent dimcnsionlcss ratios. 
The first ratio is the so-called adiabaticity parameter, 



H = -* E c h - s E( 6 * + b >i + w ° E b ^ 



(1) 



(2) 



which determines which of the two subsystems, electrons or phonons, is the 
fast or the slow one. In the adiabatic limit the motion of the particle 
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is affected by quasi-static lattice deformations (adiabatic potential surface). 
In contrast, in the anti- adiabatic limit a 3> 1, the lattice deformation is 
presumed to adjust instantaneously to the position of the carrier. The particle 
is referred to as a "light" or "heavy" polaron in the adiabatic or anti-adiabatic 
regimes [4]. 

The second ratio is the dimensionless EP coupling constant. Here 



appears in (small polaron) strong-coupling perturbation theory. Defining the 
polaron binding energy as e p — g 2 / w o = 9 2 ^o, the EP coupling can be 
parametrised alternatively as the ratio of polaron energy for an electron con- 
fined to a single site and the free electron half bandwidth 2Dt: 



In the limit of small particle density, a crossover from essentially free carri- 
ers to heavy quasiparticles is known to occur from early quantum Monte Carlo 
calculations [63], provided that two conditions, g > 1 and A > 1, are fulfilled. 
So, while the first requirement is more restrictive if a is large, i.e. in the anti- 
adiabatic case, the formation of a small polaron state will be determined by 
the second criterion in the adiabatic regime [64, 65]. 

Perhaps it is not surprising that standard perturbative techniques are 
less able to describe the Holstein system close to the large- to small-polaron 
crossover, where e p ~ 2Dt or e p ~ loq. In principle, variational approaches, 
that give correct results in the weak- and strong-coupling limits, could pro- 
vide an interpolation scheme. Most variational calculations, however, lead to a 
discontinuous transition in the wave function and the derivative of the ground- 
state energy, considered as a function of the coupling parameter. Clearly the 
analytical behaviour of an exact wave function may deviate considerably from 
that of a variational approximation [15]. With regard to the Holstein polaron 
problem the nonanalytic behaviour found for the adapted wave function in 
many variational approaches, see, e.g., [66] and references therein, is an artifact 
of the approximations, as we will demonstrate below for all dimensions [61]. 

Nevertheless, variational calculations are an indispensable tool for numer- 
ical work. In the next subsection we describe a variational exact diagonali- 
sation (VED) scheme [60] that does not suffer from the above drawback of 
(ground-state) non-analyticitics at the small-polaron transition. Above all, in 
contrast to finite-lattice ED, it yields a ground-state energy which is a vari- 
ational bound for the exact energy in the thermodynamic limit. As yet the 
VED technique is fully worked out for the single polaron and bipolaron prob- 
lem only. At finite particle densities the construction of the variational Hilbert 
space becomes delicate. On this account we will also outline some more general 
(robust) ED schemes, which can be applied for the calculation of ground-state 
and spectral properties of a larger class of strongly correlated EP systems. 



9 = g/u 



(3) 



A = e p /2Dt. 



(4) 
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2.2 Numerical techniques 

Hilbert space and basis construction 

The total Hilbert space of the Holstcin model (1) can be written as the ten- 
sorial product space of electrons and phonons, spanned by the complete basis 
set {\b) = |e) <g> \p)} with 

N N 1 

l e > = Il II (4) ni< " e |0} e and \p)=H^=(b\r^\0) p . (5) 

i=la=],i i=l V m i,P- 

Here n ia ^ e G {0, 1}, i.e. with respect to the electrons Wannier site i might 
be empty, singly or doubly occupied, whereas we have no such restriction for 
the phonon number, mj jP G {0, . . . , oo}. Consequently e = 1, . . . , D e and 
p = I,..., Dp label basic states of the electronic and phononic subspaces 
having dimensions D e = )( N N ) and D p = oo, respectively. \0) e / p 
denote the corresponding vacua. This also holds including electron-electron 
(e.g. Hubbard-type) interaction terms [67]. For Holstein-t-J-type models, act- 
ing in a projected Hilbert space without doubly occupied sites, we have 
= (jv^ ) ( N N Ne '") om y fiQ] - Since these model Hamiltonians commute with 
the total electron number operator N e = N e ,<r, where N e .a — X)t=i n i,<n 
and the z-component of the total spin S z = ^^2iLi( n i,^ — n i,i), the many- 
particle basis {\b}} can be constructed for fixed N e and S z . 

Variational approach 

Let us first describe an efficient variational exact diagonalisation (VED) 
method to solve the Holstein model in the single-particle subspace. For gener- 
alisation of this method to the case of two particles (bipolaron) see Ref. [68] . 

A variational Hilbert space is constructed beginning with an initial root 
state, taken to be an electron at the origin with no phonon excitations, and 
acting repeatedly with the hopping (t) and EP coupling (g) terms of the 
Holstcin Hamiltonian (see Fig. 1). States in generation L are those obtained 
by acting L times with these "off-diagonal" terms. Only one copy of each state 
is retained. Importantly all translations of these states on an infinite lattice 
are included. A translation moves the electron and all phonons j sites to the 
right. Then, according to Bloch's theorem, each eigenstate can be written as 
t/j = e'^ OL, where ol is a set of complex amplitudes related to the states in 
the unit cell j, e.g. L = 7 for the small variational space shown in Fig. 1. For 
each momentum k the resulting numerical problem is then to diagonalise a 
Hermitian Lx L matrix. The total number of states per unit cell (N st ) after L 
generations increases exponentially as (D + 1) L (note that the bipolaron has 
the same exponential dependence with only a larger prefactor). Most notably 
the error in the ground-state energy Eg decreases exponcntionally, because 
states are added in a fairly efficient order. Thus in most cases 10 4 - 10 6 basis 
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Fig. 1. Variational Hilbert space construction for the ID polaron problem. Basis 
states are represented by dots, off-diagonal matrix elements by lines. Vertical bonds 
create or destroy phonons with frequency ujo. Horizontal bonds correspond to elec- 
tron hops (oc t). Accordingly, state |1) describes an electron at the origin (0) and no 
phonon, state |2) is an electron and one phonon both at site 0, |3) is an electron at 
the nearest-neighbour site 1, and a phonon at site 0, and so on [60]. 



states are sufficient to obtain an 8-16 digit accuracy for Eq (see Fig. 2). The 
ground-state energy calculated this way is variational for the infinite system. 




Fig. 2. Fractional error A in the ground-state energy of a D-dimensional polaron 
as a function of the number of basis states N a t retained. Parameters are A = 0.5, 
g = 1, and t = 1. Figure is taken from Ref. [61]. 



Symmetrisation and phonon truncation 

Treating more complex many-particle Hamilton operators on finite lattices, 
the dimension of the total Hilbert space can also be reduced. To this end we 
can exploit the space group symmetries [translations (Gt) and point group 
operations (Gl)] and the spin rotational invariance [(Gg); S z — subspace 
only]. Working, e.g., on finite ID or 2D bipartite clusters with periodic bound- 
ary conditions (PBC), we do not have all the symmetry properties of the un- 
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derlying ID or 2D (square) lattices [39]. Restricting ourselves to the ID non- 
equivalent irreducible representations of the group G{K) = Gt x Gl{K) x G$, 
we can use the projection operator Vk,ts = X)seG(ff) Xk]ts & (with 

[H,V K ,rs] = 0, V\ C rs = V K ,rs and V K ,rs V K ',r' s ' = V K , rs Sk,k> &ry &s,s') 
in order to generate a new symmetrised basis set: {\b)} — > {\b)}. Q denotes the 
g(K) elements of the group G{K) and Xk\s * s tne (complex) character of Q 
in the [K, rs] -representation, where K refers to one of the N allowed wave 
vectors in the first Brillouin zone, r labels the irreducible representations of 
the little group of K, Gl(K), and s parameterises Gs- For an efficient parallel 
implementation of the matrix vector multiplication (see below) it is extremely 
important that the symmetrised basis can be constructed preserving the ten- 
sor product structure of the Hilbcrt space, i.e., 

{\b)=Nf rs] V K ,rs [|e)®|p>]} (6) 

with e = l,...,D 9 e {K) [D 9 e (K) - D e /g(K)}. The are normalisation 

factors. 

Since the Hilbcrt space associated to the phonons is infinite even for a 
finite system, we use a truncation procedure [38] retaining only basis states 
with at most M phonons: 

N 

{\p) ; m p = ^m iiP <M}. (7) 
»=i 

Then the resulting Hilbert space has a total dimension D = D^ K ^ x D^ 1 
with = (M + N)\/M\N\, and a general state of the Holstcin model is 
represented as 

\1>K,r.)= E iZ C t\b)- (8) 

e=l p=l 

The computational requirements can be further reduced if one separates the 
symmetric phonon mode, B a = -^J^i^, and calculates its contribution to 
H analytically [69]. 

Note that switching from a real space representation to a momentum space 
description the truncation scheme takes into account all dynamical phonon 
modes, which has to be contrasted with the frequently used single-mode ap- 
proach [70] . In other words, depending on the model parameters and the band 
filling, the system "decides" by itself how the M phonons will be distributed 
among the independent Einstein oscillators related to the TV Wannier sites 
or, alternatively, among the different phonon modes in Q-space. Hence with 
the same accuracy phonon dynamical effects on lattice distortions being quasi- 
localised in real space (such as polarons, Frenkel excitons,. . . ) or in momentum 
space (like charge-density-waves,. . . ) can be studied. 
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Of course, one has carefully to check for the convergence of the above 
truncation procedure by calculating the ground-state energy as a function of 
the cut-off parameter M. In the numerical work below convergence is assumed 
to be achieved if E is determined with a relative error less than 1CP 6 . 

Phonon basis optimisation 

In this section we outline an advanced phonon optimisation procedure based 
on controlled density-matrix basis truncation [56]. The method provides a 
natural way to dress the particles with phonons which allows the use of a very 
small optimal basis without significant loss of accuracy. 
Starting with an arbitrary normalised quantum state, 

IV>> = E E c&[|e>®|p>], (9) 

e=0 p=0 

expressed in terms of the basis of the direct product space, we wish to reduce 
the dimension D p of the phonon space H p by introducing a new basis, 

D p -1 

\p) = E OrpIp) . ( 10 ) 

with p — 0, . . . , (Dp — 1) and Dp < D p . We call {\p}} an optimised basis, if 
the projection of \tp) on the corresponding subspace H = H e (g> Hp C H is as 
close as possible to the original state. Therefore we minimise 

\M)-\j>)f = \-T:{ a pct) (11) 

with respect to the ap p under the condition (p'\p) = 5p' P , where 

D e -1 Dg-1 D p -1 

w = E E E ^p^4e[\^)®\p)\ (12) 

e=0 p=0 p,p'=0 

is the projected state, p = Y^=o^ ( c t P Y ' c t P * s called the density matrix of 
the state \ip). Clearly the states {\p}} are optimal if they are elements of the 
eigenspace of p corresponding to its Dp largest eigenvalues w p . If we interpret 
w p ~ exp(— ap) as the probability of the system to occupy the corresponding 
optimised state \p) , we immediately find that the probability for the complete 
phonon basis state ^^S^ 1 \pi) is proportional to exp(— a X^o* Pi)- This is rem- 
iniscent of an energy cut-off, and we therefore propose the following choice of 
a mixed phonon basis at each site, 

Vz: {| Mi >} = ON({| M )}) (13) 
Joptimal state \p), < p < M opt 
~~ { bare state \p), M opt < p, < M ' ' 
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sweep 1 sweep 2 

Fig. 3. Sweep technique in constructing optimised phonon states. 



and for the complete phonon basis {®£ <w <Af|Mi)}) yielding D p h = ( N+ ^~ 1 )- 
After a first initialisation the optimised states are improved iteratively 
through the following steps: (1) calculate the requested eigenstate of the 
Hamiltonian H in terms of the actual basis, (2) replace {\p}} with the most 
important (i.e. largest eigenvalues wp) eigenstates of the density matrix p; (3) 
change the additional states {\p}} in the set (4) orthonormalise the set 

{|£t)}, and return to step (1). 

A simple way to proceed in step (3) is to sweep the bare states {\p)} 
through a sufficiently large part of the infinite dimensional phonon Hilbert 
space. One can think of the algorithm as "feeding" the optimised states with 
bare phonons, thus allowing the optimised states to become increasingly per- 
fect linear combinations of bare phonon states (see Fig. 3). 



Solution of the eigenvalue problem 



In all the above cases the numerical problem that remains is to find the eigen- 
states of a (sparse) Hermitian matrix. Here iterative (Krylov) subspace meth- 
ods like Lanczos [51] and variants of Davidson [71] diagonalisation techniques 
are frequently applied. These algorithms contain basically three steps: (1) 
project the problem matrix A S R™ onto a subspace A fc S V fc (i; < n); (2) 
solve the eigenvalue problem in V fc using standard routines; (3) extend the 
subspace V fe — ► V fc+1 by a vector t 1 V 1 and go back to (2). This way we 
obtain a sequence of approximate inverses of the original matrix A. A pow- 
erful and widely used technique is the Lanczos algorithm which recursively 
generates a set of orthogonal states (Lanczos vectors): 



\tpt +1 )=E D \(pi)-ai\<pi)-^\ipi- 1 ) ) (15) 



where 



_ (<pi\H. D \<pi) ,2 _ (<Pi\<Pi) ,2_ n nfn 

ai - -, : r , 0; - -. : r , - U , (lb) 

WW) W-iW-i) 

and \ip~i) = 0. Obviously, the representation matrix [T L ]/j/ = ((pi\H D \tpi>) 
of H D is tridiagonal in the L-dimensional Hilbert space spanned by the 
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{|<^;)}/=o,...,L-i ; where L <C D. The eigenvalue spectrum of T L can be eas- 
ily determined using standard routines from libraries such as EISPACK (see 
http://www.netlib.org). Note that the convergence of the Lanczos algorithm 
is excellent at the edges of the spectrum (the ground state for example is ob- 
tained with high precision after at most ~ 100 Lanczos iterations) but rapidly 
worsens inside the spectrum. 

In general the computational requirements of these eigenvalue algorithms 
are determined by matrix- vector multiplications (MVM), which have to be 
implemented in a parallel, fast and memory saving way on modern supercom- 
puters. The MVM step can be be done in parallel by using a parallel library 
such as PETSc (see http://www-unix.mcs.anl.gov/petsc/petsc-as/). 

Our matrices are extremely sparse because the number of non-zero entries 
per row of our Hamilton matrix scales linearly with the number of electrons. 
Therefore a standard implementation of the MVM step uses a sparse storage 
format for the matrix, holding the non-zero elements only. The typical storage 
requirement per non-zero entry is 12-16 Byte, i.e. for a matrix dimension of 
D = 10 9 about one TByte main memory is required to store only the matrix 
elements of the EP Hamiltonian. To extend our EP studies to even larger ma- 
trix sizes we no longer store the non-zero matrix elements but generate them 
in each MVM step. Of course, at that point standard libraries are no longer 
useful and a parallel code tailored to each specific class of Hamiltonians must 
be developed. Clearly the parallelisation approach follows the inherent natural 
parallelism of the Hilbert space. Assuming that the electronic dimension (D e ) 
is a multiple of the number of processors used we can easily distribute the 
electronic basis states among the processors. As a consequence of this choice 
only the electronic hopping term generates inter-processor communication in 
the MVM while all other (diagonal electronic) contributions can be computed 
locally on each processor. Using supercomputers with hundreds of processors 
and one TByte of main memory, such as the IBM p690 cluster, at the moment, 
one is able to run simulations up to a matrix dimension of 3 x 10 10 . 

Determination of dynamical correlation functions 

The numerical calculation of dynamical spectral functions, 



A°(lo) = - lim -Im 

D-1 



= J2 l(^|0|Vo)| 2 ^ - (E n - Eo)} , (17) 



n=0 



where O is the matrix representation of a certain operator O [e.g., the 
creation operator <J K of an electron with wavevector K if one wants to 
calculate the single-particle spectral function; or the current operator j = 
— iet X^( c ! c i+i — c !+i c i) ^ onc i s interested in the optical conductivity], in- 
volves the resolvent of the Hamilton matrix H. Once we have obtained the 
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eigenvalues and eigenvectors of H we can plug them into Eq. (17) and directly 
obtain the corresponding dynamical correlation or Green functions. For the 
typical EP problems under investigation we deal with Hilbert spaces having 
total dimensions D of 10 6 - 10 11 . Finding all eigenvectors and eigenstates of 
such huge Hamilton matrices is impossible, because the CPU time required 
for exact diagonalisation of H scales as D 3 and memory as D 2 . So in practice 
this "naive" approach is applicable for small Hilbert spaces only, where the 
complete diagonalisation of the Hamilton matrix is feasible. Fortunately, there 
exist very accurate and well-conditioned linear scaling algorithms for a direct 
approximate calculation of A°(ui). 

Kernel polynomial method (KPM) 

The idea of the KPM (for a review see [53]), is to expand A°(u>) in a finite 
series of L + 1 Chebyshcv polynomials T m (x) = cos[m arccos(x)]. Since the 
Chebyshev polynomials are defined on the real interval [—1,1], first a sim- 
ple linear transformation to the Hamiltonian and all energy scales has to be 
applied: X = (H - b)/a, x = (w - b)/a, a = (E max - E min )/2(1 - e), and 
b = (E max + E m i n )/2 (the small constant e is introduced in order to avoid 
convergence problems at the endpoints of the interval - a typical choice is 
e ~ 0.01 which has only 1% impact on the energy resolution [72]). Then the 
expansion is 

A°{x) = J__ L° + 2 jr &T m (x)\ , (18) 

X \ m=l / 

with moments 

»Z = j 1 ^dxT m (x)A°{x) = {^T m {X)0\^). (19) 

Eq. (18) converges to the correct function for L — > oo. The moments 

V% m = 2 (4>m\^m) - Mo and $ m+l = 2((j) m+1 \4> m ) - nf (20) 

can be efficiently obtained by repeated parallelised MVM, where \4> m +i) = 
2X\(f> m ) — \4> m -i) but now \<px) — X|O0o) with \(j> ) = \tpo) determined by 
Lanczos ED. 

As is well known from Fourier expansion, the series (18) with L finite 
suffers from rapid oscillations (Gibbs phenomenon) leading to a poor ap- 
proximation to A°(oj). To improve the approximation the moments fi n are 
modified fi n — ► g n ^n, where the damping factors g n arc chosen to give the 
"best" approximation for a given L. In more abstract terms this truncation 
of the infinite series to order L together with the corresponding modification 
of the coefficients is equivalent to a convolution of the spectral function with 
a smooth approximation kernel Ki,(x,y). The appropriate (optimal) choice 
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of this kernel, that is of g n , e.g. to guarantee positivity of A°(uj), lies at the 
heart of KPM. We mainly use the Jackson kernel which results in a uniform 
approximation whose resolution increases as 1/L, but for the determination 
of the single-particle Green functions below we use a Lorentz kernel which 
mimics a finite imaginary part r\ in Eq. (17) [53]. 

In view of the uniform convergence of the expansion, the accuracy of the 
spectral functions can be made as good as required by just increasing L. 

Cluster perturbation theory ( CPT) 

The spectrum of a finite system of TV sites which we obtain through KPM 
differs in many respects from that in the thermodynamic limit N — > oo, espe- 
cially in that it is obtained for a finite number of momenta K = itm/N only. 
While we cannot easily increase N without going beyond computationally ac- 
cessible Hilbert spaces we can try to extrapolate from a finite to the infinite 
system. 

For this purpose we first calculate the Green function G? (w) for all sites 
i, j = 1, . . . , N c of a N c - size cluster with open boundary conditions, and then 
recover the infinite lattice by pasting identical copies of this cluster at their 
edges [54]. The "glue" is the hopping V between these clusters, where Vki = t 
for \k — l\ = 1 and k, I = 0, l(modiV), which is dealt with in first order pertur- 
bation theory. Then the Green function Gijiui) of the infinite lattice is given 
through a Dyson equation Gij{w) = G?-(w) + J2ki Gi k (u)VkiGij(uj), where 
indices of G c (lo) are counted modulo N c . The Dyson equation is solved by 
Fourier transformation over momenta K — kN c corresponding to translations 
by N c sites 



G ij {K,u) = 
from which one finally obtains 



V(K)G c {w) 



(21) 



N c 

G(k,cj) = G^.(Ayc, W )exp(-ifc(i- j)) . (22) 

In this way we obtain a Green function G(k, lu) with continuous momenta 
k from the cluster Green function G c . Two approximations are made, one 
by using first order perturbation theory in V = t, the second on assuming 
translational symmetry in Gij (u>) which is only approximatively satisfied. 



3 Ground state results 

The VED method can compute polaron properties in dimensions 1 through 

4 and higher. The energies for ID to 4D polarons at k = for intermediate 
to weak coupling on linear, square, cubic, and hypercubic lattices are listed 
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Table 1. Polaron ground state energies at k — in ID - 4D for A = 0.5 and g = 1.0 
In the following t = 1 unless specifically noted. 



ID 


2D 


3D 


4D 


E /t -2.46968472393287071561 


-4.814735778337 


-7.1623948409 


-9.513174069 



in the Tab. 1. The number of significant digits is determined by comparing 
the energy as the size of the Hilbert space is increased. The accuracy is high 
compared to that of other numerical methods, even when limited by the single- 
processor desktop workstations of several years ago, on which the results were 
obtained [61]. Correlation functions are generally less accurate than energies. 




Fig. 4. Logarithm of the polaron effective mass in ID as a function of g. VED results 
(full lines) were obtained operating repeatedly L — 20 times with the off-diagonal 
pieces of the Holstein Hamiltonian (cf. Sec. 2.2). For comparison GL (global-local) 
results (dashed lines) are included [47]. Open symbols, indicating the value of coo/t, 
are DMRG results [48]. 



Figure 4 shows the effective mass computed by VED [60] in comparison 
with GL and DMRG methods, m* is obtained from 

mo = 1 d 2 E k 
to* ~~ It dk 2 



(23) 



Numerical solution of the Holstein polaron problem 



15 



where m = l/{2t) is the effective mass of a free electron. The second deriva- 
tive is evaluated by small finite differences in the neighbourhood of k = 0. 
Note that although the calculated energy E k is a variational bound for the 
exact energy, there is no such control on the mass, which may be either above 
or below the exact answer, and is expected to be more difficult to obtain 
accurately. Nevertheless, in the intermediate coupling regime where VED at 
L = 20 gives an energy accuracy of 12 decimal places, one can calculate the 
effective mass accurately (6-8 decimal places) by letting Ak — > 0. 

In Fig. 4 the parameters span different physical regimes including weak and 
strong coupling, and low and high phonon frequency. We find good agreement 
between VED and GL away from strong coupling and excellent agreement in 
all regimes with DMRG results. DMRG calculations are not based on finite-fc 
calculations due to a lack of periodic boundary conditions, so they extrap- 
olate the effective mass from the ground state data using chains of different 
sizes, which leads to larger error bars and demands more computational effort. 
Notice that their discrete data is slightly scattered around the VED curves. 
Nevertheless, both methods agree well. We have compared the VED results for 
effective mass obtained on different systems from L = 16 with N st = 178617 
states to L = 20 with N st — 2975104 states and obtained convergence of re- 
sults to at least 4 decimal places in all parameter regimes presented in Fig. 4. 
Our error is therefore well below the lincwidth. Even though there is no phase 
transition in the ground state of the model, the polaron becomes extremely 
heavy in the strong coupling regime. The crossover to a regime of large polaron 
mass is more rapid in adiabatic regime, i.e. at small u)o/t. 

The polaron effective mass in higher dimensions is shown as a function of 
the EP coupling in Fig. 5. The mass increases exponentially for large A. The 
crossover to larger effective mass is more rapid, though still continuous, in 
higher dimensions. 

Of course it is of special interest to understand the evolution of the pola- 
ronic band structure as the EP coupling strength increases. Figure 6 presents 
the results for the wavevector dependence of the ground-state energy Ek in ID 
at weak (a) and strong (b) EP couplings. As might be expected, for A = 0.25 
the coherent bandwidth, AE = sup fc Ek — inffc Ek, is approximately given by 
the phonon energy (AE = 0.782 <~ uj = 0.8). If the EP interaction is en- 
hanced a band collapse appears. Note, however, that even in the relatively 
strong EP coupling regime displayed Fig. 6 (b) the standard Lang-Firsov for- 
mula, AElf = 4Dt cxp[— g 2 } (obtained by performing the Lang-Firsov canon- 
ical transformation [30] and taking the expectation value of the kinetic energy 
over the transformed phonon vacuum), does not give a satisfactory estimate 
of the bandwidth. So we found AElf = 0.0269 which has to be contrasted 
with the exact result AE = 0.0579. Besides the band narrowing effect, there 
are several other features worth mentioning for polaronic band states in the 
crossover region. Most notably, throughout the whole Brillouin zone the band 
structure differs significantly from that of a rescaled tight-binding (cosine) 
band containing only nearest-neighbour hopping [65]. Obviously the residual 
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Fig. 5. Effective mass of the Hoistein poiaron in dimensions 1, 2, and 3. 



polaron-phonon interaction generates longer-ranged hopping processes [3, 65]. 
Concomitantly, the mass enhancement due to the EP interaction is weakened 
at the band minimum. It is important to realize that these effects are most 
pronounced at intermediate EP couplings and phonon frequencies. In this pa- 
rameter region even higher-order strong-coupling perturbation theory, with 
its internal states containing some excited phonons, seems to be intractable 
because the convergence of the series expansion is poor [73]. Of course the 




Fig. 6. Band structure of the ID Hoistein model in the weak (a) and strong (b) EP 
coupling regimes. M and L denote the number of phonons and the depth of the basis 
in ED and VED calculations, respectively. Within ED the (Q = 0) centre of mass 
motion has been separated if indicated. The inset shows the finite-size dependence 
being significant for weak EP coupling only. In the strong EP coupling case ED 
results basically agree with those obtained by VED. 
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Fig. 7. The quasiparticle weight Zk and in the inset the total number of phonons 
N£ h as a function of the wavevector k for the ID Holstein model with (g/t) 2 = 0.4 
[dashed line (VED); open circles (ED)], 3.2 [solid line (VED); filled circles (ED)]. 
a — LOo/t = 0.8. Data are taken from Refs. [60, 65]. 



dispersion is barely changed from a rescaled tight-binding band in the very 
extreme small polaron limit (A, g 2 3> 1). 

Further information about the quasiparticle may be obtained by comput- 
ing the quasiparticle residue, the overlap squared between a bare electron and 
a polaron, 

^ = K^l4lO>| 2 , (24) 

where |0) is the state with no electron and no phonon excitations, and |V>fc) 
is the polaron wavefunction at momentum k. Zj~ can be measured by angle- 
resolved photoemission, and gives the bare electronic contribution of the pola- 
ronic state. The phonon contribution to the quasiparticle is described by the 
k— dependent mean phonon number 



(25) 



Figure 7 shows the spectral weight Zp. and the mean phonon number TV? 
as a function of k. The two sets of parameters correspond to the large and 
small polaron regime respectively. The DMRG cannot straightforwardly com- 
pute this quantity. There is excellent agreement between the VED and ED 
methods in the weak coupling case. In the strong coupling regime there is 
an approximately 1% disagreement in N% h due to a lack of phonon degrees 
of freedom in the variational space of the ED calculation. The results in the 
weak coupling case show a smooth crossover from predominantly electronic 
character of the wavefunction for small k (large Zp. and small N^ h w 0) to 
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predominantly phonon character around k = n characterised by Z k w and 
N? w 1. In the strong coupling regime there is less fc-dependence. The Z^ is 
already close to zero at small k, indicating strong EP interactions that lead to 
a large polaron mass. Concomitantly an appropriately defined polaron quasi- 
particle residue Zq tends to one [74, 75]. So we arrive at the conclusion that 
a well-defined electronic (polaronic) quasiparticle exists for k = at weak 
(strong) EP coupling. 

VED is one of the few methods that can also calculate the polaron band 
dispersion in 3D systems (QMC is another, but is limited to the condition 
that the polaron bandwidth is much smaller than the phonon frequency, which 
corresponds to the strong-coupling regime.) Figure 8(a) shows the band dis- 
persion for the 3D polaron along symmetry directions in the Brillouin zone at 
various EP coupling constants g = gu . Starting with weak coupling g = 2.0 
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Fig. 8. Ground-state energy Eh (shifted by Eo) of the 3D polaron in panel (a) 
and quasiparticle weight in panel (b) for three different EP coupling constants, 
g — 4.5 (solid line), g — 3.5 (dotted line), and g = 2.0 (dashed line), for uio/t = 2. The 
dot-dashed line in (a) is the dispersion of a bare electron. The corresponding ground 
state energies E /t are -10.608348, -8.0642850, and -6.588526818 respectively. 



(dashed line), the polaron band is close to the bare electron band at the 
lower band edge. The deviation between them increases as k increases. When 
Ek — Eq approaches u)q, we observe a band flattening effect, similar to the ID 
and 2D cases, accompanied by a sharp drop of quasiparticle weight Z]~ (see 
Fig. 8(b)). The large k lowest energy state can be considered roughly as "a 
k = polaron ground state" plus "an itinerant (or weakly-bound) phonon 
with momentum fc". It is the phonon that carries the momentum so as to 
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make Zf. essentially vanish and give an approximate bandwidth ujq. Due to 
the large spatial extent of the EP correlations in the flattened band, our results 
are less accurate in this regime. In the case of intermediate coupling g = 3.5 
(A = 1.0208), the polaron bandwidth is narrower than the phonon frequency. 
The upper part of the band has much less dispersion than the lower part but 
with a substantial Zk- This indicates a distinct mechanism for the crossover 
as a function of k. In the case of g = 4.5, the strong EP interaction leads to a 
significant suppression of Zk at all k. Zk=a approaches the strong-coupling re- 
sult Z = exp(— g 2 ) for A , g 2 ^> 1. Note that the inverse effective mass m* /mo 
and Z differ if the self-energy is strongly fc-dependent. This discrepancy has 
its maximum in the intermediate-coupling regime for ID systems, but van- 
ishes in the limit A — > oo. In the Holstein model with on-site electron- phonon 
interactions, Z and the inverse effective mass are closely related. However, 
the two can be made arbitrarily different by increasing the range of the EP 
interaction [61]. 

The correlation function between the electron position and the phonon 
displacement is 

X id = (^ k \clc t (b 3 +b])\Tjj k ). (26) 

This correlation function can be considered as a measure of the polaron size. 
It should not be confused with the "polaron radius" in the extreme adiabatic 
limit, which refers to the spatial extent of a hypothetical symmetry-breaking 
localised state. The ground-state EP correlation function is plotted for the 
ID Holstein polaron in Fig. 9. For parameters close to the (adiabatic) weak 
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Fig. 9. Renormalised EP correlation function xo,j = ( n o(&o +:? + b 0+j )) /2g{no) as a 
function of the electron-phonon separation j in the k = ground state of the ID 
Holstein model. Results are taken from [17, 76]. 



EP coupling regime (main panel), the amplitude of Xo,j is small and the 
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spatial extent of the electron-induced lattice deformation is spread over the 
whole (finite) lattice. The DMRG results shown in the left inset indicate a 
substantial reduction of the polaron's size near the crossover region. Finally 
a small polaron is formed at large couplings (right inset); now the position of 
the electron and the phonon displacement is strongly correlated. 

How does the electron-phonon correlation function change as the polaron 
acquires a nonzero momentum kl The answer is shown for ID in Fig. 10 [60]. 
The parameters correspond to weak coupling. At k = 0, where the group 
velocity is zero, the deformation is limited to only a few lattice sites around 
the electron. The correlation is always positive and exponentially decaying. 
At finite but small k = 7r/4, the deformation around the electron increases in 
amplitude and rings (oscillates in sign) as the polaron acquires a finite group 
velocity. At k = ir/2 the ringing is strongly enhanced. Note also that the 
spatial extent of the deformation increases in comparison to k — 0. The range 
of the deformation is maximum at k — 7r, where it extends over the entire 
region shown in the figure. In keeping with the larger extent of the lattice 
deformation near k = 7r, the ground-state energy E„ converges more slowly 
with the size of the Hilbert space. We have also computed the k— dependent \ 
for the strong-coupling case uj = 0.8, g 2 = 3.2 (not shown). We find only weak 
fc— dependence, which is a consequence of the crossover to the small polaron 
regime. The lattice deformation is localised predominantly on the electron 
site. 

Over four decades ago, a simple and intuitive variational approach to the 
ID polaron problem was proposed by Toyozawa [28]. This method has been 
successfully applied to various fields and revisited in a number of guises over 
the years. It is generally believed to provide a qualitatively correct description 
of the polaron ground state, aside from predicting a spurious discontinuous 
change in the mass at intermediate coupling. The Toyozawa variational wave- 
function consists of a product of coherent states (displaced oscillators) around 
the instantaneous electron position. The phonons create a symmetrical cloud 
around the electron. Numerical studies of the ID electron-phonon correlation 
function (two-point function) are in semi-quantitative agreement with the 
Toyozawa variational wavefunction. The numerically exact three-point func- 
tion, however, disagrees wildly. Denoting the instantaneous electron position 
as 0, the Toyozwa variational wavefunction requires that the probability to 
find phonon excitations, for example, on both sites 3 and 4, should be iden- 
tical to finding them on sites (-3) and 4. Numerically, however, the latter is 
many orders of magnitude less probable [61]. This suggests a physical picture 
in which the k = polaron is viewed not as an electron surrounded by a 
symmetrical cloud of phonons, but is instead a coherent superposition of two 
"comets," one with a tail extending to the right, and the other to the left. 

Studying the properties of Holstein polarons, DMFT is exact in infinite di- 
mensions. An interpolation to 3D lattices is made possible by using a semiellip- 
tical free density of states N(E) to mimic the low-energy features. Here DMFT 
is accurate in the strong-coupling regime, where the surrounding phonons are 
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predominately on the electron site. This is also the regime where strong- 
coupling perturbation theory works well. However, DMFT fails to compute 
quantities such as the polaron mass correctly in the weak-coupling regime. The 
reason is that in DMFT, the lattice problem is mapped onto a self-consistent 
local impurity model, which preserves the interplay of the electron and the 
phonons only at the local level. The spatial extent of the EP correlations 
increases as the EP coupling decreases, which explains the significant discrep- 
ancy in the weak-coupling regime. Therefore only the on-site EP correlation 
has been studied by DMFT, and the results are compared with VED [61] in 
Fig. 11. There is good qualitative agreement. The curves show a rapid change 
in slope only for large g, where DMFT is less accurate. It is worth noting that 
DMFT neglects the k dependence of self-energy, i.e., the inverse effective mass 
is always equal to the quasiparticle spectral weight. Clearly the difference be- 
tween m / m* and the spectral weight Zk is not negligible in the intermediate- 
to weak-coupling regime. 
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4 Excited states 



In this section we turn our attention from the ground state to the excited 
states of the Holstein model. Figure 12 plots the energy eigenvalues for a 
small variational space containing a maximum of 9 phonon excitations. The 
lowest curve is the polaron ground state at momentum k. Excited states are 
the polaron with additional bound or unbound (or both) phonon excitations. 
A ripple can be discerned near the bare electron dispersion. The figure su- 
perficially resembles a "band structure," which however encodes ground and 
excited state information for the many-body (many phonon) polaron problem. 
The ac conductivity of the polaron, for example, appears as an "interband" 
transition in this mapping. 

What is the nature of the first excited state? We focus on the question of 
whether an extra phonon excitation forms a bound state with the polaron, or 
instead remain as two widely separated entities. Using numerical and analyt- 
ical approaches we show evidence that there is a sharp phase transition (not 
a crossover) between these two types of states, with a bound excited state 
appearing only as the EP coupling is increased. Although the ground state 
energy E is an analytic function of the parameters in the Hamiltonian, there 
are points at which the energy E\ of the first excited state is nonanalytic. 
In previous work, Gogolin has found bound states of the polaron and addi- 
tional phonon(s), but he does not obtain a phase transition between bound 
and unbound states because his approximations are limited to strong coupling 




Fig. 11. On-site correlation x(0) f° r the 3D polaron. VED results (solid lines) [61] 
are compared to DMFT (dotted lines with symbols) [34]. x(0) is normalised to 1 
when A is infinite (i.e. t — > 0) according to lim x(0) = 2g. Parameters are uo = t = 1. 



Numerical solution of the Holstein polaron problem 



23 




0.0 0.2 0.4 0.6 0.8 1.0 

k 



Fig. 12. The ground and excited state ID polaron energy eigenvalues (those energies 
below 0) are plotted as a function of k (in units of 7r) for g = u>o = t = 1, L = 9, 
Nst = 1185. 




Fig. 13. First excited state binding energy A = E\ — Eq — uio as a function of g [60] . 
Results are for lo — 0.5 and various Hilbert space sizes L. Inset: Binding energy 
over a wider range of g. 

3 > 1 [77]. A phase transition between a bound and unbound first excited 
state has been calculated for 3D using a dynamical CPA (coherent potential 
approximation) [33] and DMFT [34]. 
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We compute the energy difference AE = E x — E , where E x and E are 
the first excited state and the ground state energies at k — (the two lowest 
bands in Fig. 12). In the case where the first excited state of a polaron can be 
described as a polaron ground state and an unbound extra phonon excitation, 
this energy difference should in the thermodynamic limit equal the phonon 
frequency, AE = uq. In Fig. 13 we plot the binding energy A = AE — ojo for 
ojq = 0.5 as a function of the EP coupling g for various sizes of the variational 
space. We see two distinct regimes. Below g c ~ 0.95, A varies with the system 
size but remains positive (A > 0). Physically, for g < g c , the additional 
phonon excitation would prefer to be infinitely separated from the polaron, 
but is confined to a distance no greater than L — 1 by the variational Hilbcrt 
space. As the system size increases, A slowly approaches zero from above 
as the "particle in a box" confinement energy decreases. In the other regime, 
g > g c , the data has clearly converged and A < 0. This is the regime where the 
extra phonon excitation is absorbed by the polaron forming an excited polaron 
bound state. Since the excited polaron forms an exponentially decaying bound 
state, the method already converges at L = 14. In the inset of Fig. 13 we 
show the binding energy A in a larger interval of EP coupling g. Although 
the results cease to converge at larger g, we notice that the binding energy 
A reaches a minimum as a function of g. As one can demonstrate within 
the strong coupling approximation, the true binding energy approaches zero 
exponentially from below with increasing g. 

Figure 14 shows the phase diagram for k = separating the two regimes. 
The phase boundary, given by A = 0, was obtained numerically, and compared 
to strong coupling perturbation theory in t to first and second order. The phase 
transition where A becomes negative at sufficiently large g is also seen in ED 
calculations. 

The distribution of the number of phonons in the vicinity of the electron 
is defined as 

7(i-j) = ^fc|ctc i 6t6 J .|V fc >. (27) 

In Fig. 15 we compute this distribution for the ground state 70 and the first 
excited state 71 slightly below (g = 0.9), and above (g — 1.0) the transition 
for luq = 0.5. 

The central peak of the correlation function 71 below the transition point 
is comparable in magnitude to 70 (Fig. 15 (a,b)). The main difference between 
the two curves is the long range decay of 71 as a function of distance from the 
electron, onto which the central peak is superimposed. The extra phonon that 
is represented by this long-range tail extends throughout the whole system and 
is not bound to the polaron. The existence of an unbound, free phonon is con- 
firmed by computing the difference of total phonon number Nq x = ^ 7o,i(0- 
This difference should equal one below the transition point. Our numerical 
values give Nf h — N^ h <~ 1.02. We attribute the deviation from the exact re- 
sult to the finite relative separation allowed. Correlation functions above the 
transition point (Fig. 15 (c,d)) are physically different. First, phonon corre- 
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Fig. 14. The phase diagram for the bound to unbound transition of the first excited 
state, obtained using the condition A(u>o,g) = 0. The corresponding phase diagram 
for the ground state would be blank: there is no phase transition in the ground state, 
only a crossover. 

lations in 71 decay exponentially, which also explains why the convergence in 
this region is excellent. Second, the size of the central peak in 71 is 3 times 
higher than 70. (Note that to match scales in Fig. 15 (d) we divided 71 by 
3). The difference in total phonon number gives 

N ph _ N ph _ 2 33 We 

are 

thus facing a totally different physical picture: The excited state is composed 
of a polaron which contains several extra phonon excitations (in comparison 
to the ground-state polaron) and the binding energy of the excited polaron 
is A < 0. The extra phonon excitations are located almost entirely on the 
electron site. The value of 71 — 70 at j = is 2.16, which almost exhausts the 
phonon sum. 

Next we discuss the role of dimensionality in the excited states. The effect 
of dimensionality on static properties has been studied previously [78, 79, 80, 
61]. The eigenvalues of the low- lying k — states are shown as functions of 
g in Fig. 16. The energy spectra in D>1 arc qualitatively different than in 
ID. The ID polaron ground state becomes heavy gradually as g increases. 
However, in D>2, the ground state crosses over to a heavy polaron state by 
a narrow avoided level crossing, which is consistent with the existence of a 
potential barrier [78]. In the lower panel of Fig. 16, -01 and 04 are nearly 
free electron states; -02 and -03 are heavy polaron states. The inner product 
KV'il'i/^)! is equal to 0.99. Just right of the crossing region the effective mass 
(approximately equal to the inverse of the spectral weight) of the first excited 
state can be smaller than the ground state by 2 or 3 orders of magnitude, while 
their energies can differ by much less than uj . The narrow avoided crossing 
description works less well for larger uj . 
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Fig. 15. The phonon number 7 as a function of the distance from the electron 
position for the ground state (a) and the first excited state (b), both computed at 
g = 0.9; and the same in (c) and (d) for g = 1.0. All data are computed at phonon 
frequency loo = 0.5 and L — 18. Note that (d) is a plot of 71/3. In (a,b), 71 — 70 
drops to zero around \i — j\ = 15. This is a finite-size effect. Computing the same 
quantity with larger L below the phase transition would lead to a larger extent of 
the correlation function indicating that the extra phonon excitation is not bound to 
the polaron. 



5 Dynamics of polaron formation 

How does a bare electron time evolve to become a polaron quasiparticle? The 
bare electron can be injected by inverse photocmission or tunnelling, or a hole 
by photocmission, or an exciton (electron-hole bound state) by fast optics. 

One approach is to construct a variational many-body Hilbcrt space in- 
cluding multiple phonon excitations, and to numerically integrate the many- 
body Schrodingcr equation, 

i d f t =m (28) 

in this space [81]. The full many-body wave function is obtained. This method 
includes the full quantum dynamics of the electrons and phonons. Note that 
alternative treatments, such as the semiclassical approximation that treats the 
phonons classically, fail for this problem, particularly in the limit of a wide 
initial electron wavepacket. 

Figure 17 shows snapshots of polaron formation at weak coupling. An 
initial bare electron wave packet is launched to the right as shown in panel 
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Fig. 16. Eigenvalues of low-lying states as functions of coupling constant in ID 
through 3D. Hopping t = 1 in all panels. In the adiabatic regime in higher dimen- 
sions, the ground state (thick solid lines) shows a fairly abrupt change in slope. In 
the 3D panel, ipi and ip4 are a lightly-dressed electron state; V>2 and ip3 are a heavy 
polaron state. The dashed lines are the beginning of the lowest continua. 



(a). This initial condition is relevant to the recent experiments [82, 83, 84, 85, 
86], and to electron injection from a time-resolved STM (scanning tunneling 
microscope) tip [87]. Although polarons injected optically or by STM can 
have a range of initial momenta, it would be more realistic to take k = 
for an optically created exciton. In panel (b) the electron is not yet dressed 
and thus is moving roughly as fast as the free electron (green dashed line). 
In addition, there exists a back-scattered current (which later evolves into a 
left-moving polaron) on the left side of the wave packet (green dot-dashed and 
thick black curves). In panel (c) after an elapsed time of one phonon period, 
the electron density consists of two peaks. The peak on the right (black arrow) 
is essentially a bare electron. The peak on the left is a polaron wave packet 
moving more slowly. As time goes on, the bare electron peak decays and the 
polaron peak grows. Some phonons are left behind (blue double-dot dashed 
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line), mainly near the injection point. These phonons are of known phase with 
displacement shown in thin solid red. Some phonon excitations travel with the 
polaron (magenta dot double-dashed line). Finally, a coherent polaron wave 
packet is observed when the polaron separates from the uncorrected phonon 
excitations. The velocity operator is defined as 

.-,. = 2 km (29) 



e(c] Cj + c] +1 c j+i y 

where j is the site index and j is the current operator, (ij) is shown as a green 
dot-dashed line. 

There are regimes where the polaron formation time is a calculable con- 
stant of order unity times a phonon period T , as seen in some experiments 
and in Fig. 17, but there are other regimes in which the phonon period is 
not the relevant timescale. The limit of hopping t — > is instructive [88, 89]. 
After a time To/4, the expectation of the lattice displacement (xj) on the 
electron site has the same value as a static polaron. It is tempting (but we 
would argue incorrect) to identify this as the polaron formation time. At later 
times, (xj) overshoots by a factor of two, and after a time To, (xj) and all 
other correlations are what they were at time zero when the bare electron was 
injected. The system oscillates forever. In general an electron emits phonons 
enroute to becoming a polaron, and we propose that the polaron formation 
time be defined as the time required for the polaron to physically separate 
from the radiated phonons. The polaron formation time for hopping t — > 
is thus infinite, because the electron is forever stuck on the same site as the 
radiated phonons. 

An electron injected at several times the phonon energy ujq above the 
bottom of the band is another instructive example. The electron radiates 
successive phonons to reduce its kinetic energy to near the bottom of the 
band, and then forms a polaron. For very weak EP coupling, the rate for 
radiating the first phonon can be computed by Fermi's golden rule, TpQ R — 
g 2 /[fr isin(fc/)], where kf is the electron momentum after emitting a phonon. 
The phonon emission time can be arbitrarily long for small g. For strong 
coupling, the rate approaches t<7 c = g/K because the polaron spectral function 
smoothly spans numerous narrow bands and its standard deviation is equal 
to g. 

Decaying oscillations in polaron formation (actually the formally equiva- 
lent problem of an exciton coupled to phonons [4]) have been observed in a 
pump-probe experiment that measures reflectivity after a bare exciton is cre- 
ated [83]. The observed oscillatory reflectivity was interpreted as the lattice 
motion in the phonon-dressed (or "self-trapped") exciton level. Assuming the 
modulation in the exciton level goes as AE = —XcjCjXj, where Xj = (bj + &t) 
is the lattice displacement, the model Hamiltonian applies directly to the ex- 
periment. We calculate the corresponding EP correlation function in Fig. 18. 
In this regime, the polaron formation time (damping time) increases as the 
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Fig. 17. Snapshots of the polaron-formation process, for t = Uq = 1, and g — 0.4. 
The calculation is performed on a 30-site periodic lattice. Time is measured in units 
of the phonon period. Shown are the electron density (c^Cj) (thick solid black line), 
phonon density (blue double-dot dashed line), lattice displacement (xj) = 

(bj +bj) (thin solid red line), velocity in units of lattice constant per phonon period 
(green dot dashed line), and EP correlation function (c^Cjbjbj) (magenta dot double- 
dashed line). The green dashed line gives the free-electron wave packet for reference. 
For clarity, the origins of the thin solid red and green dashed curves are offset by 
0.1 and their values are rescaled by a factor of 0.2 and 0.05/(27r) respectively. The 
blue double-dot dashed curve has been rescaled by a factor of 0.5. 

electron-phonon coupling g increases, and also as the initial electron (exciton) 
energy approaches the band bottom. We find satisfactory agreement when 
compared to Fig. 2b of Ref. [83]. Both show a damped oscillation with a 
delayed phase. (Numerical calculations in Figs. 18-19 are performed on an 
extended system, not a finite cluster.) 

Figure 19 shows the spectral function at strong coupling. Three delta func- 
tions are visible, corresponding to polaron ground and excited states, along 
with three continua containing unbound phonons. There is additional struc- 
ture at higher energy (not shown). The probability to remain in the initial 
bare particle state P(t) = |(V , ( r )l c I|0)| 2 for this spectrum is complicated, 




Fig. 18. The on-site electron-phonon correlation function \ — {c^Cj{bj + b^)) as a 
function of time measured in phonon periods. For all curves, luq — 0.5 and hopping 
t = 1. The solid line is for a bare electron injected with nonzero initial momentum at 
energy Ei = —0.7, where the bottom of the bare band is at energy -2. The phonon 
displacement is larger and more weakly damped for larger electron-phonon coupling 
g, dotted line. In contrast to a bare electron, an exciton (bound particle-hole pair) is 
generally created with an initial momentum zero, corresponding to Ei = —2, dashed 
line. 



and includes oscillating terms that do not decay to zero at zero temperature 
from the polaron ground and excited states beating against each other. The 
branching ratios into the various channels are calculated in Ref. [90]. 



6 Spectral signatures of Holstein polarons 

As already stressed in the introduction the crossover from quasi-free elec- 
trons or large polarons to small polarons becomes manifest in the spectral 
properties above all. Here of particular interest is whether an "electronic" or 
"polaronic" (quasi-particle) excitation exists in the spectrum. This question 
has been partially addressed by calculating the wavefunction renormalisation 
factor [(electronic) quasi-particle weight] Zk in Sec. 3 (see Fig. 7). More de- 
tailed information can be obtained from the one-particle spectral function 
A{k,w). This quantity of great importance can be probed by direct (inverse) 
photoemission, where a bare electron with momentum k and energy lo is re- 
moved (added) from (to) the many-particle system. The intensities (transition 
amplitudes) of these processes are determined by the imaginary part of the 
retarded one-particle Green's function, i.e. by 




Fig. 19. Panel (a): Spectral function at strong coupling. There are three quasipar- 
ticle excited states split off from the continua. Shaded areas (1) and (2) correspond 
to continuum states, (b) : Quantum beat formed by multiple excited states and con- 
tinua. 



A{k,ui) = --hnG(k,w) = A+(k,Lu) + A~ (k,w) , (30) 

7T 



with 



A ± (k,uj) = Imlim (^ |c^ — — — cj|^. 



7T rj->0+ oj + 17] =F H 

= £ T (E± - E )] , (31) 

m 

2^ = Cj. and CjT — c k (T = 0; ID spinless case). These functions test 
both excitation energies (£"„ — E ) and overlap (c>c Zk) of the A^-particle 
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ground state \tp ) with the exact eigenstates |^) of a (N e ± l)-particle sys- 
tem. The electron spectral function of the single-particle Holstein model cor- 
responds to N e — 0, i.e., A~(k,uj) = 0. A(k,u) can be determined, e.g., by a 
combination of KPM and CPT (cf. Sec. 2.2). 

Figure 20 (a) shows that at weak EP coupling, the electronic spectrum is 
nearly unaffected for energies below the phonon emission threshold. Hence, for 
the case considered with ujo lying inside the bare electron band = — 2t cos k, 
the signal corresponding to the renormalised dispersion Ek nearly coincides 
with the tight-binding cosine band (shifted oc e p ) up to some kx, where the 
phonon intersects the bare electron band. At kx electron and phonon states 
"hybridise" and repel each other, forming an avoided-crossing like gap. For k > 
kx the lowest absorption signature in each k sector follows the dispersionless 
phonon mode, leading to the flattening effect [65]. Accordingly the (electronic) 
spectral weight of these peaks is very low. The high-energy incoherent part 
of the spectrum is broadened oc e p , with the fc-dependent maximum again 
following the bare cosine dispersion. 

Reaching the intermediate EP coupling (polaron crossover) regime a co- 
herent band separates from the rest of the spectrum [kx — > tt; see panel (b)]. 
At the same time its spectral weight becomes smaller and will be transferred 
to the incoherent part, where several sub-bands emerge. 

The inverse photoemission spectrum in the strong-coupling case is shown 
in Fig. 20 (c). First, we observe all signatures of the famous polaronic band- 
collapse. The coherent quasi-particle absorption band becomes extremely nar- 
row. Its bandwidth approaches the strong-coupling result 4Diexp(— g 2 ) for 
A,g 2 S> 1. If we had calculated the polaronic instead of the electronic spec- 
tral function, nearly all spectral weight would reside in the coherent part of 
the spectrum, i.e. in the small-polaron band. This has been demonstrated 
quite recently [75]. In our case the incoherent part of the spectrum carries 
most of the spectral weight. It basically consists of a sequence of sub-bands 
separated in energy by loq, which correspond to excitations of an electron and 
one or more phonons. 

Let us emphasise that for all couplings the lowest signature in A(k, to) 
almost perfectly coincides with the coherent polaron band-structure (solid 
lines) obtained by VED (sec Sec. 3), which underlines the high precision of 
the CPT-KPM approach used here. 

Of course, the phonon modes are unaffected by one electron in the solid, i.e. 
the phonon self-energy is zero. Actually this is true in the thermodynamic limit 
only. In a finite-cluster calculation there will be an influence of order 0(1/N) 
and the phonon spectra provide additional useful information concerning the 
polaron dynamics. For this purpose, we calculate the T = phonon spectral 
function 

B(q,u) = ImD R (q,uj) (32) 

which is related to the phonon Green's function 
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Fig. 20. Spectral function of the ID Holstein polaron in the weak (a), intermediate 
(b), and strong (c) EP coupling regimes. CPT data based on finite-cluster ED with 
N c sites, and M = 7 (A = 0.25), M = 15 (A = 1), M = 25 (A = 2) phonon quanta. 
Note that the non-monotonic heights of the lowest energy peaks in (a) are an artifact 
of the CPT calculation, where some of the wavevectors fit the N c = 16 cluster size, 
and some don't. Also the dispersionless absorption feature in (c), just above the 
small-polaron peak, is due to a finite-size effect, but not the double-peak structures 
of the higher excitation bands. This has been proved by determining the spectral 
function in the k = 0-sector by VED. 
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D R (q,uj)= lira (i/j \x q — — — x_ 5 |i/>o) , (33) 

»)^o+ u; + it] — H 

where x q = N' 1 ' 2 £\ ^e _i « and x 3 = (b\ + 

For the Holstein model (1), B{q 1 Lj) is symmetric in q. The bare propaga- 
tor D (q,ui) = 2ui / (ui 2 — ujq) is dispersionless. Then, adapting the CPT-KPM 
approach to the calculation of the phonon spectral function, the cluster ex- 
pansion is identical to replacing the full real-space Green's function Dij by 
D c . 




Fig. 21. Electron (a) and phonon (b) spectral functions in the anti-adiabatic inter- 
mediate EP coupling regime. Solid (dashed) lines give Ek (u>o) determined by VED. 
Note that abscissae are scaled differently. 



Figure 21 compares electron (a) and phonon (b) spectra in the high 
phonon- frequency limit, where the small polaron crossover is determined by 
g 2 . Obviously the phonon spectrum is also renormalised by the EP interaction 
due to the finite "particle density" N e /N c = 1/10. So we can detect a clear 
signature of the polaron band having a width W ~ 1.5t (cf. Fig. 21 (a)). The 
dispersionless excitation at lu/loq — 1 is obtained by adding one phonon with 
momentum q to the k = ground state. Above this pronounced peak, we find 
an "image" of the lowest polaron band - shifted by Wo ~ with extremely small 
spectral weight, hardly resolved in Fig. 21 (b). 



7 Transport and optical response 

The investigation of transport properties has been playing a central role in 
condensed matter physics for a long time. Optical measurements, for example, 
proved the importance of EP interaction in various novel materials such as the 
cuprates, nickelates or manganites and, in particular, corroborated polaronic 
scenarios for modelling their electronic transport properties at least at high 
temperatures [91, 92, 93]. Actually, the optical absorption of small polarons is 
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distinguished from that of large (or quasi-free) polarons by the shape and the 
temperature dependence of the absorption bands which arise from exciting 
the self-trapped carrier from or within the potential well that binds it [16]. 
Furthermore, as was the case with the properties of the ground state, the 
optical spectra of light and heavy electrons, small and large polarons differ 
significantly as well [4]. In the most simple weak-coupling and anti-adiabatic 
strong EP coupling limits, the absorption associated with photoionization 
of Holstein polarons is well understood and the optical conductivity can be 
analysed analytically ([91, 94, 95, 96, 66]; for a detailed discussion of small 
polaron transport phenomena we refer to Refs. [97, 98]). The intermediate 
coupling and frequency regime, however, is as yet practically inaccessible for a 
rigorous analysis (here the case of infinite spatial dimensions, where dynamical 
mean-field theory yields reliable results, is an exception [99, 100]). So far 
unbiased numerical studies of the optical absorption in the Holstein model 
were limited to very small 2 to 10-site ID and 2D clusters [36, 64, 74, 17]. In 
the following we will exploit the VED and KPM schemes [62, 101], in order to 
calculate the optical conductivity numerically in the whole parameter range 
on fairly large systems. 



7.1 Optical conductivity at zero-temperature 

Applying standard linear-response theory, the real part of the conductivity 
takes the form 

RccrH = VS(u) + a reg (w) , (34) 

where V denotes the so-called Drude weight at u> — and a reg is the regular 
part (finite-frequency response) for uo > which can be written in spectral 
representation at T = as [94] 

are9{uj) = ^N £ \(^m\Mo)\ 2 S[u - (E m - E Q )} (35) 

with the (paramagnetic) current operator j = — iet X)i( c l c i+i — c \+i c i)- 
Introducing the w-integrated spectral weight, 

S re9 (cu) = f du'a re9 (u'), (36) 
Jo+ 

we arrive at the f-sum rule 



-E kin /2 = V + S reg /ir (ID case), 



(37) 



where Eu n = —t Si( c l c i+i+ c l+i c i) * s t ne kinetic energy and S reg — S reg (oo). 
Since for the Holstein model the Drude weight can be calculated independently 
from Kohn's formula or the effective mass, 



V 



1 d 2 E (<P) 



2N d$ 2 



1 d 2 E k 



4>=a 2N dk 2 



k=0 



2m* 



(38) 
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the f-sum rule may be used to test the numerics. In Eq. (38), the first equality 
relates V to the second derivative of the (non-degenerate) ground-state energy 
with respect to a field-induced phase coupled to the hopping. 

We first present a reg (uj) and its integral S reg {uo) for the ID Holstein model 
in Fig. 22. The upper panel (a) gives the results for intermediate-to-strong 
EP coupling, i.e. near the polaron crossover, in the adiabatic (light electron) 
regime. Of course, the optical conductivity threshold is u = ujo for the infinite 
system we deal with using VED. In this respect standard ED, defined on finite 
lattices, suffers from pronounced finite-size effects due to the discreteness in k- 
space. Knowing that at about A ~ 1 a coherent polaron band with bandwidth 
much smaller than luq splits off, the first (few) isolated peak(s) at the lower 
bound of the spectrum can be attributed to one- (few-) (q = 0— ) phonon 
emission processes (cf. also Fig. 20 (b)). Of course, these transitions have 
little spectral intensity. At higher energies transitions to the incoherent part 
of the spectrum take place by "emitting" phonons with finite momentum (to 
reach the total momentum k = ground-state sector). The main signature of 
a re9 {iij) is that the spectrum is strongly asymmetric, which is characteristic 
for rather large polarons. The weaker decay at the high-energy side meets the 
experimental findings for many polaronic materials like Ti02 [102] even better 
than standard small-polaron theory. 

For A = 2 and uo = 0.4, i.e., at larger EP coupling, but not yet in the 
small-polaron limit, we find a more pronounced and symmetric maximum in 
the low-temperature optical response (see Fig. 22 (b)). The maximum is lo- 
cated below the corresponding one for small polarons at T = 0, which on 
its part lies somewhat below 2e p = 2t\ = 2g 2 u; (being the maximum of the 
Poisson distribution) because of the 1/oj factor contained in the conductiv- 
ity. In this case the polaron band structure is more strongly renormalised, 
but, more importantly, the phonon distribution function in the ground state 
becomes considerably broadened. Since the current operator connects only 
different-parity states having substantial overlap as far as the phononic part is 
concerned, in the optical response multi-phonon emissions/absorptions (i.e., 
non-diagonal transitions [94]) become increasingly important. Again devia- 
tions from the analytical small-polaron result (dashed line in Fig. 22 (b)) 
might be important for relating theory to experiment. 

The optical response obtained if the phonon frequency becomes compara- 
ble to the electron transfer amplitude is illustrated in Fig. 22 (c). Now the 
lowest one-phonon absorption (threshold) signal is well separated. In contrast 
to the light electron case Fig. 22 (a), the different absorption bands appearing 
for a heavier electron can be classified according to the number of phonons 
involved in the optical transition (see inset). Increasing w at fixed g 2 this 
becomes even more manifest (at the same time a Poisson distribution of the 
different sub-bands evolves). The sub-bands are broadened by transitions to 
different "electronic" levels. For our parameters, a scattering continuum ap- 
pears above ui > 5 to 6 uq. Note that the "fragmentation" of the spectrum 
appearing at smaller energy transfer is not caused by finite-size effects. 
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Fig. 22. Optical conductivity of theoJP2S°l s tein polaron at T = 0. Results for 

a reg {u) and S reg {u) are obtained by VED. I?i order to reduce finite-M effects, data 

calculated for M = 15 to 20 were averaged. The dashed line in (b) displays the 



analytical small polaron result, a reg (ui) 



"a 1 



exp 



(cf. Ref. [91]), 



where ctq was determined to give the same integrated spectral weight as a re9 (uj > 0). 
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Turning to the sum rules presented in Fig 23, we notice a monotonic de- 
crease of the total sum rule S tot /ir = —E^ in /2, which indicates a suppression 
of the electronic kinetic energy with increasing EP coupling. In agreement 
with previous numerical results [63, 17, 43], the kinetic energy clearly shows 
the crossover from a large polaron, characterised by a E^in that is only weakly 
reduced from its non-interacting value [Eki n {^ = 0) = — 2t], to a less mobile 
small polaron in the strong EP-coupling limit, where the strong-coupling per- 
turbation theory result, 

fcm LU \s/ K=2g2 

((. . .) K denotes the average with respect to the Poisson distribution with pa- 
rameter k), gives a sufficiently accurate description in both the adiabatic and 
antiadiabatic regimes. 



4i / 1 



LOq \ S I K=g 



(39) 




Fig. 23. Ronormalised kinetic energy (Ekin, solid line) and contribution of a reg to 
the f-sum rule (S reg ; dot-dashed line) as a function of the EP couplings A and g in the 
adiabatic (left panels) and non-to-antiadiabatic (right panels) regimes, respectively. 
The Drude weights were obtained from the f-sum rule [(Eq. (37); thin solid line)] 
and effective mass [(Eq. (38); dashed line)]. 



For light electrons (adiabatic regime io^jt — 0.2, 0.4; left panels), we found 
a rather narrow transition region. The drop of S tot in the crossover region 
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A ~ 1 is driven by the sharp fall of the Drude weight, which is a measure of the 
coherent transport properties of a polaron. By contrast the optical absorption 
due to inelastic scattering processes, described by the regular (dissipative) part 
of the optical conductivity, becomes strongly enhanced around A ~ 1 [74] (cf. 
the behaviour of S re9 ). 

The large to small polaron crossover is considerably broadened for heavy 
electrons (non-to-antiadiabatic case u>o/t = 1, 4; right panels). Here Ekin 
decreases more gradually and S reg exhibits a less pronounced maximum at 
about g 2 = 1. 

As quoted above, we can calculate the Drude weight independently from 
the effective mass of the Holstein polaron. Using this data, it is worth men- 
tioning that the f-sum rule (37) is satisfied numerically to at least six digits 
in the whole parameter regime [62] . 

7.2 Thermally activated transport 

If the polaron effects are assumed to be dominant the coherent bandwidth is 
extremely small. Then the physical picture is that the particle is trapped at 
a certain lattice site and that hopping occurs infrequently from site to site. 
There are two kinds of transfer processes [11]. All phonon numbers might 
remain the same during the hop (diagonal transition) or, alternatively, the 
number of phonons is changed (non-diagonal transition). In the latter case 
each hop may be approximated as a statistically independent event and the 
particle loses its phase coherence by this phonon emission or absorption (in- 
elastic scattering). Diagonal and non-diagonal transitions show a different 
temperature dependence. While the rate of diagonal (band-type) transitions 
decreases with increasing temperature, small-polaron theory predicts that the 
non-diagonal (incoherent hopping) rate is thermally activated and may be- 
come the main transport process at higher temperatures (cf., e.g., Rcf. [94]). 
Deviations from standard small-polaron theory are expected to occur in the 
intermediate coupling regime. By means of ED and KPM techniques we are 
able to study the optical response of Holstein polarons precisely in this regime, 
at least for small lattices. 

ac conductivity 

Our starting point is the Kubo formula for the electrical conductivity at finite 
temperatures [94], 

oo 

ff ™>;T) = JL- £ [ e -^--e-^»] \W> n \jW m )\ 2 6{u>-w mn ), (40) 

m,n>0 

where Z = e~ /3B ™ is the partition function and (3 = T~ x denotes the 
inverse temperature (ftg = 1)- Since in practice the contribution of highly ex- 
cited phonon states is negligible at the temperatures of relevance, the system 
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is well approximated by a truncated phonon space with at most M(A, g, uj a ; T) 
phonons [39]. Then \ip n ) and \ip m ) are the eigenstates of H within our trun- 
cated Hilbert space. E n and E m are the corresponding eigenvalues with 

<^mn — E m E n . 

In order to evaluate temperature-dependent response functions like (40), 
recently a generalised "two-dimensional" KPM scheme has been proposed [53, 
101], which, in our case, can be set up using a current operator density 

j(x,y)=Y,\^\^\ 2 S ( x ~ E n) 5{y-E m ). (41) 

m.n 

For the regular part of the conductivity we obtain 

oo 

° rC9 ^) = ^\ Jj(y + u,y)[e-ev-e-^ + ^]dy, (42) 
o+ 

where the partition function Z = 2 J °^ p(E) exp(— /3.E) is easily obtained by 

integrating over the density of states p(E) = J2n=o ~ E n ), which can 
be expanded in parallel to j(x,y) (here D is the dimension of the Hilbert 
space). One advantage of this approach is that the current operator density 
that enters the conductivity is the same for all temperatures, i.e., it needs to 
be expanded only once. 

Figure 24 gives results for the finite-temperature optical conductivity of 
a Holstcin polaron. Coherent transport related to diagonal transitions within 
the lowest polaron band is negligible at high temperatures. For instance, the 
amplitude of the current matrix elements between the degenerate states with 
momentum K = ±7r/3 (K — 0, ±7r/3, ±27r/3, and 7r are the allowed wave 
numbers of our 6-site system with periodic boundary conditions) is of the 
order of 10~ 7 only. In the small polaron limit , where the polaronic sub-bands 
are roughly separated by the bare phonon frequency, non-diagonal transitions 
become important for T > uj . Let us consider the activated regime in more 
detail (see Fig. 24 (upper panel)). With increasing temperatures we observe a 
substantial spectral weight transfer to lower frequencies, and an increase of the 
zero-energy transition probability in accordance with previous results [103]. 

In addition, we find a strong resonance in the absorption spectra at about 
ujq ~ 2i, which can be easily understood using a configurational coordinate 
picture [101]. In order to activate these transitions thermally, the electron has 
to overcome the "adiabatic" barrier A = E\ + — Eq = e p /2 — t, where we have 
assumed that the first relevant excitation is a state with lattice distortion 
spread over two neighbouring sites and the particle mainly located at both 
these sites (in a symmetric (+) or antisymmetric (— ) linear combination; 
E\^± = =Fi — £ p /2). A finite phonon frequency will relax this condition. From 
Fig. 24, we find the signature to occur above T > 0M. Obviously this feature 
is absent in the standard small-polaron transport description which essentially 
treats the polaron as a quasiparticle without resolving its internal structure. 
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Fig. 24. Optical absorption by Holstein polarons at finite temperatures in the 
adiabatic strong (upper panel) and intermediate (lower panel) EP coupling regime. 
Results are obtained by ED for a TV = 6 site lattice with M = 45 phonons. In the 
upper panel, thin lines with symbols give the analytical results for the small polaron 
transport [104, 94] at temperatures f3t = 2 (triangles), 0.5 (squares), 0.2 (circles). 
The deviations observed for high excitation energies at very large temperatures are 
caused by the necessary truncation of the phonon Hilbert space in ED. 
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Now let us decrease the EP coupling strength A keeping g 2 = 10 fixed. 
Results for the optical response in the vicinity of the large to small polaron 
crossover are depicted in the lower panel of Fig. 24. Here the small polaron 
maximum has almost disappeared and the 2t-absorption feature can be acti- 
vated at very low temperatures (A — > for the two-site model with A = 1). 
The gap observed at low frequencies and temperatures is clearly a finite-size 
effect. The overall behaviour of a re9 (u}; T) resembles that of polarons of inter- 
mediate size. At high temperatures these polarons will dissociate readily and 
the transport properties are equivalent to those of electrons scattered by ther- 
mal phonons. Let us emphasise that many-polaron effects become increasingly 
important in the large-to-small polaron transition region [105] (see also Sec. 8 
below). As a result, polaron transport might be changed entirely compared to 
the one-particle picture discussed so far. 

dc conductivity and thermopower 

We consider dc transport, or a(uj) in the limit w — > 0. For simplicity, we con- 
sider only a single polaron, or a dilute system of polarons where interactions 
can be neglected and bipolaron formation is prevented, as by a large repulsive 
U. We also neglect impurities, which can localise or scatter a polaron. 

At zero temperature, the conductivity or mobility of a polaron is infinite. 
The polaron can be placed in a state of nonzero momentum by a weak electric 
field acting for a short time. This is an eigenstate, which carries current forever 
and never decays. At small temperatures T < wo, an exponentially small 
number of phonons are thermally excited. The conductivity becomes finite 
due to scattering of a polaron off thermally excited phonons of density n p h ~ 
e -^o/T _ rp ne details depend on the EP scattering process. 

In ID, when a polaron of momentum k encounters a thermally excited 
phonon, in general part of it is transmitted and part is backscattered. Certain 
anomalies occur. For example, in the limit of small hopping t, as g approaches 
1, the backscattering of the polaron vanishes and the phonon is simultaneously 
transferred one site in the direction opposite the polaron momentum. The 
phonon thus recoils opposite to the direction expected, cf. the collision of two 
balls. This leads to a heat current in the opposite direction as the polaron 
particle current, which should be observable in the thermopower. A polaron- 
thermal phonon bound state also exists for sufficiently large g. For this bound 
state, heat (a phonon excitation) can be transported by an electric field, which 
again should be observable as a large contribution to the thermopower of the 
opposite sign as the above. For large g, this bound state or internal polaron 
excited state can have a much smaller effective mass than the polaron ground 
state. Perhaps surprisingly, as the temperature increases, the polaron effective 
mass as measured by the low-frequency ac conductivity can decrease. 

We next consider very high temperatures. As T increases, the typical 
phonon displacement increases as K(x 2 ) = ksT, where K is the phonon 
spring constant. For quasi-static phonons (large phonon mass), this leads to a 
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disorder potential for the electron that increases without bound as T increases. 
The disorder Anderson localises the electron, leading to zero dc conductivity. 
The disorder, however, is not quite static, and rearranges itself on a timescale 
t ~ 1/wo- Once every time of order t, the diagonal energies of the electron 
site and a neighbouring site become equal, and the electron can hop to a 
neighbouring site. It is then diffusing with a diffusion constant <~ a 2 uja, where 
a is the lattice constant. Using the Einstein relation relating diffusion and 
mobility, the high temperature resistivity becomes 

P=^T2 — • ( 43 ) 

The high temperature resistivity is metallic, i.e. dp/dT > 0, and can greatly 
exceed the Ioffe-Regel limit. Numerical studies to confirm or refute this sce- 
nario are incomplete. 



8 From few to many polarons 

Let us now address the important issue of how the character of the (pola- 
ronic) quasiparticles may change if we increase the carrier density n = N e /N. 
Consider first the case of zero electron-electron interaction. Beginning with 
a noninteracting Fermi gas at T = 0, as the Holstein EP interaction g is 
increased from zero, a singlet superconductor is expected to form. As g in- 
creases, the diameter of the Cooper pair decreases. Eventually, the Cooper 
pair diameter becomes smaller than the distance between Cooper pairs, and 
the behaviour crosses over from BCS superconductivity to that of Bosc con- 
densation, like that of 4 He, where the hard core bosons are bipolarons (bound 
states of two polarons). In this limit, T c is given approximately by the Bose 
condensation temperature for ideal bosons of mass m* , where m* is the bipo- 
laron mass. The limit of Bose condensation of bipolarons is not given correctly 
by Eliashberg theory, which describes strong coupling, but not that strong. 

8.1 Bipolaron formation 

We investigate how two electrons coupled to phonons may bind together to 
form a bipolaron, including the bipolaron effective mass, the crossover between 
two different types of bound states, and the dissociation into two polarons 
(see also [106, 107]). For problems with more than one electron, the Holstein 
Hamiltonian is generalised by adding a Hubbard electron-electron interaction 
term, U^2j n^n^. Basis states for the many-body Hilbert space can be writ- 
ten \b) = \ji,j2'i ■ ■ ■ ,n m ,n m+ i, ...,), where the up and down electrons are on 
sites ji and j2, and there are n m phonons on site m. In a generalisation of 
the one electron VED method described above, a bipolaron variational space 
is constructed beginning with an initial state where both electrons are on the 
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same site with no phonons, and operating repeatedly (L-times) with the off- 
diagonal pieces (t and g) of the Hamiltonian. All translations of these states 
are included on an infinite lattice. The method is very efficient in the inter- 
mediate coupling regime, where it provides results that are variational in the 
thermodynamic limit and bipolaron energies that are accurate to 7 digits for 
the case L = 18 and size of the Hilbert space N st — 2.2 x 10 6 phonon and 
down electron configurations for a given up electron position. In ID the size 
of the variational space approximately doubles as L is increased by one, which 
is the same as for the one electron problem, although the prefactor for two 
electrons is larger [68]. 

For large phonon frequency u>o, the EP interaction leads to a non-retarded 
attractive on-site interaction of strength J7o = 2uj g 2 - One would expect that 
as the Hubbard repulsion U becomes larger than this value, the bipolaron 
would dissociate into two polarons. As can be shown both analytically and 
numerically, this is not what happens. In the limit of small hopping t, as U 
exceeds Uo, the bipolaron crosses over from a state SO with both electrons 
primarily on the same site, to another bound state SI with the electrons pri- 
marily on nearest neighbour sites. Only for U > 2Uo does the bipolaron disso- 
ciate into two polarons. The crossover from SO to 51 bipolarons is important 
in theories of bipolaronic superconductivity applied to real materials, since 
SI bipolarons arc generally orders of magnitude lighter than SO bipolarons. 
Since the superconducting T c in the dilute limit is inversely proportional to 
the effective mass, the SI regime usually provides a more compelling theory. 

We now discuss numerical variational results for the singlet bipolaron on 
an infinite ID lattice. We have been unable to demonstrate the existence 
of a bound triplet bipolaron for the Holstein-Hubbard model. Fig. 25 shows 
the ground state electron-electron density correlation function C(i — j) = 
(ipo\n>iiij\ipo) , where n, = mi + mi and |-0o) denotes the ground state wave 
function. At g = 1, the bipolaron widens with increasing U and transforms 
into two unbound polarons (which can only move a finite distance apart in 
the variational space). The value U — 1.5 is below the transition to the un- 
bound state at U c = 2.17 ', calculated by comparing the polaron and bipolaron 
energies. We see that the probability of electrons occupying the same or neigh- 
bouring sites is almost equal. In the unbound regime, the nature of the corre- 
lation function changes significantly. At U = 1.5, C(i) falls off exponentially, 
while for U > U c the typical distance between electrons is the order of the 
maximum allowed separation L. The electrons can be no farther apart than 
L in the variational space, although their centre of mass can be anywhere on 
an infinite lattice. A state of separated polarons is clearly seen for U — 20. 

Two distinct regimes are seen at g = 2 within the bipolaronic region. At 
U = 7 < Uq = 2a>o.9 2 = 8) the correlation function represents the SO bipolaron, 
while at U — 9 > Uo we find the largest probability for two electrons to be 
on neighbouring sites, which is characteristic of the SI bipolaron. In contrast 
to previous calculations where phonons were treated classically [108], we find 
a crossover rather than a phase transition between the two regimes. Moving 
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from strong towards intermediate coupling, the SO and SI bipolarons consist 
of longer exponential tails extending over many lattice sites, and the two 
regimes can no longer be distinguished. The precision of presented correlation 
functions in the bipolaron regime is within the size of the plot symbols in the 
thermodynamic limit. 
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Fig. 25. Electron-electron correlation function C(i — j) calculated at coo = 1, a) 
g = 1 and b) g = 2 for different values of U, with L = 18. The two ordinate axes 
have a different range. Insets show results for U = 0. All curves are normalised, 
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Figure 26a plots the bipolaron mass ratio R m — mbi/2m po vs. U for differ- 
ent values of uj and g. In all cases presented in Fig. 26, R m approaches 1 as 
U approaches U = U c in agreement with a state of two free polarons. At fixed 
uiq = 1 the bipolaron mass ratio increases by several orders of magnitude with 
increasing g at U = 0. Increasing U sharply decreases Rm in the SO regime. 
Note that the mass scale is logarithmic. In the SI regime with U > Uo, R m is 
small, as predicted by the strong coupling result. 

In the dilute bipolaron regime, the bipolaron isotope effect is the same 
as the classic superconductivity isotope effect for T c . The bipolaron isotope 
effect, shown in Fig. 26b, is large in the strong coupling (w = 1,5 = 2) and 
small U regime, where its value is somewhat below the large g strong coupling 
prediction a/, so ~ 2g 2 — j = 7.75. With increasing U, a/ decreases and in the 
SI regime approaches a/ ; si = g 2 /2 = 2. A kink is observed in the crossover 
regime. With decreasing g or uj 0} a.i also decreases. 
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Fig. 26. a) The mass ratio R m = mbi/2m po vs. U and b) the bipolaron isotope 
effect ai vs. U. Numerical results are for L = 18. Results for R m at loq = 0.5 
are obtained by extrapolating L — > oo. Precision in all curves is within the line- 
width in the thermodynamic limit, except for ai with uio = 0.5, where the error is 
estimated to be ±5%. The thin line in (a) is the strong coupling expansion result 
for uio = 1, g = 2. Polaron masses in units of the noninteracting electron mass are 
m po = 1.35, 1.76, 10.4, 3.06 from top to bottom respectively. 



The phase diagram U c (g) is shown in Fig. 27 at fixed too — 1- Numerical 
results, shown as circles, indicate the phase boundary between two dissociated 
polarons each having energy E po and a bipolaron bound state with energy ■ 
In the inset of Fig. 27 we show the bipolaron binding energy defined as A = 
Em — 2E po . The phase diagram is obtained from A = 0. The dashed line, given 
by Uq — 2cjo.9 2 j is a reasonable estimate for the phase boundary at small g. 
At large g the dashed line roughly represents the crossover between a massive 
SO and lighter SI bipolaron. The SI region grows with increasing g. The dot- 
dashed line is the phase boundary between SI and the unbound polaron phase, 
as obtained by degenerate strong coupling perturbation theory. Numerical 
results approach this line at larger g. The dot-dashed line asymptotically 
approaches U c — 4w <7 2 - 

8.2 Many-polaron problem 

We consider how polarons evolve from the dilute to the concentrated limit, in 
the regime where spinless or fully spin-polarised fermions prevent bipolaron 
formation. In ID with open boundary conditions, the spinless fermion prob- 
lem is equivalent to infinite Hubbard U. While for very strong EP coupling no 
significant changes are expected due to the existence of rather independent 
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Fig. 27. Phase diagram and binding energy A in units of t (inset) calculated at 
loo = I- Numerical results are circles. For greater accuracy, results near the weak 
and strong coupling regime were obtained by extrapolating L — > co. 



small (self-trapped) polarons with negligible residual interaction (assuming 
spinlcss fcrmions or strong enough electron-electron repulsion to prevent bipo- 
laron formation), a density-driven crossover from a state with large polarons 
to a metal with weakly dressed electrons should occur in the intermediate- 
coupling regime. This issue has recently been investigated theoretically by 
ED [109], QMC [105], and variational canonical transformation [75] methods, 
and is known to be of experimental relevance, e.g., in La 2 /3(Sr/Ca)i/ 3 Mn03 
films [110]. 

In the spinless fermion Holstein model, the above-mentioned density- 
driven transition from large polarons to weakly EP-dressed electrons is ex- 
pected to be possible only in ID, where large polarons exist at intermediate 
coupling. The situation is different for Frohlich-type models [7, 111, 112] with 
long-range EP interaction, in which large-polaron states exist even for strong 
coupling and in D > 1. 

To set the stage, we first comment on the evolution of the one-electron 
spectral function A(k,u>) with increasing electron density n in the weak- and 
strong-coupling limiting cases [105]. In the former the spectra bear a close 
resemblance to the free-electron case for all n, i.e., there is a strongly disper- 
sive band running from — 2t to 2t, which can be attributed to weakly dressed 
electrons with an effective mass close to the non-interacting value. As ex- 
pected, the height (width) of the peaks increases (decreases) significantly in 
the vicinity of the Fermi momentum. In the opposite strong-coupling limit the 
spectrum exhibits an almost dispersionless coherent polaron band V n < 0.5. 
Besides, there are two incoherent features located above and below the Fermi 
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energy, broadened oc e p , which are due to phonon- mediated transitions to 
high-energy electron states. The most important point, however, is the clear 
separation of the coherent band from the incoherent parts even at large n, 
indicating that small polarons are well-defined quasiparticles in the strong- 
coupling regime, even at high carrier density. 

Figure 28 displays the inverse photoemission L4 + (fc, u)] and photoemission 
spectra [A~(k, u>)] at intermediate EP coupling strength, determined by CPT. 
At low densities, n = 0.1 (upper panel), we can easily identify a (coherent) 
polaron band crossing the Fermi energy level Ep = fi(T — > 0), the latter 
being situated at the point where A~ and A + intersect. This large-polaron 
band has rather small electronic spectral weight especially away from Ep and 
flattens at large k, as known from single-polaron studies (see Sec. 6, Fig. 20). 
Below this band, there exist equally spaced phonon satellites, reflecting the 
Poisson distribution of phonons in the ground state. Above Ep there is a broad 
dispersive incoherent feature whose maximum closely follows the dispersion 
relation of free particles. 

As the density n increases, a well-separated coherent polaron band can no 
longer be identified. At about n ~ 0.3 the deformation clouds of the (large) 
polarons start to overlap leading to a mutual (dynamical) interaction between 
the particles. Increasing the carrier density further, the polaronic quasipar- 
ticles dissociate, stripping their phonon cloud. This is the case shown in the 
lower panel of Fig. 28. Now diffusive scattering of electrons and phonons seems 
to be the dominant interaction mechanism. As a result both the phonon peaks 
in A~(k,u) and the incoherent part of A + (k,cu) are washed out, the spectra 
broaden and ultimately merge into a single wide band. Most notably, the in- 
coherent excitations now lie arbitrarily close to the Fermi level. Obviously the 
low-energy physics of the system can no longer be described by single-particle 
small-polaron theory. 

9 Polaronic effects in strongly correlated systems 

The interplay of electron-electron and electron-phonon interactions in the 
formation of dressed quasiparticles is becoming the focus of attention in 
many contexts, including conducting polymers, ferroelectrics, halide-bridgcd 
transition- metal chain complexes, and several important classes of perovskitcs. 
Especially research on high-T c superconductivity (HTSC) and colossal mag- 
netoresistance (CMR) has spurred intense investigations of the competition 
or, if possible, of the cooperation of these two fundamental interactions (for 
a recent review see [26] , and references therein) . 

Many experiments have indicated substantial EP interaction in the high- 
Tc cuprates. The relevance of EP coupling can be seen from the experimen- 
tal observation of phonon renormalisation [113]. Ion channelling [114, 115], 
neutron scattering [116] and photo-induced absorption measurements [117] 
proved the existence of large anharmonic lattice fluctuations, which may be 
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Fig. 28. Single particle spectral functions A (k,u) (red dashed lines) and A + (k,u>) 
(black solid lines) for two characteristic band fillings, n = 0.1 and n — 0.4, at 
w /t = 0.4 and A = 1 (T = 0). Results are obtained by CPT using N c = 10). Blue 
crosses track the small-polaron band determined by ED. 
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responsible for local phonon-driven charge instabilities in the planar Cu0 2 
electron system [118, 119]. Photo-induced absorption experiments [120], in- 
frared spectroscopy [121] and reflectivity measurements [122] indicate the for- 
mation of small polarons in the insulating parent compounds La2Cu04 +y and 
Nd2Cu04_ y of the hole- and electron-doped superconductors La2- x Sr x Cu04+ y 
and Nd2- x Ce x Cu04_ y , respectively. Recently angle-resolved photoemission 
spectroscopy data were interpreted in terms of strong EP coupling giving rise 
to self- localisation of holes (hole polarons) [123]. Based on these experimental 
findings, several theoretical groups [124, 125, 126, 127, 128, 129, 130, 131, 132, 
133, 134] promote a (bi)polaronic scenario for HTSC. 

Even stronger evidence for polaron formation in doped charge transfer 
oxides is provided by experiments on the nickelates La2- x Sr x Ni04 [135, 23]. 
The isostructural compounds La2CuC>4 and La2NiC>4 show a remarkable dif- 
ference upon the substitution of La by Sr. Both materials become metallic 
upon doping, but in the nickelates a nearly total substitution of La for Sr is 
necessary. Also in La 2 - x Sr x Ni04 +y no superconductivity is found for any x. A 
resolution of this problem might be given by extended LDA calculations [136], 
which show that the nickelates are much more susceptible to a breathing po- 
laron instability than the cuprates. The reason is the much stronger magnetic 
confinement effect of additional holes and nickel spins. These low -spin com- 
posite holes are nearly entirely prelocalised and the EP coupling becomes 
much more effective in forming polarons. For the composition Lai.sSro.sNiC^ 
(quarter filling, x = 0.5), electron diffraction measurements show a commen- 
surate superstructure spot at the (n, 7r)-point, which has been interpreted as 
a sign of truly 2D ordering of breathing- type polarons, i.e., as a polaronic 
superlattice. 

Localised lattice distortions are also suggested to play an important role in 
determining the electronic and magnetic properties of hole-doped manganese 
oxides of the form Lai_ x [Sr, Ca] x Mn03 [137, 24]. In the region xmi ~ 0.2 < 
x < 0.5 these compounds show a transition from a metallic ferromagnetic 
low-temperature phase to an insulating paramagnetic high-temperature phase 
associated with a spectacularly large negative magnetoresistive response to 
an applied magnetic field [138]. Both breathing-mode collapsed (Mn 4+ ) and 
(anti) Jahn- Teller distorted (Mn 3+ ) sites are created simultaneously when 
the holes are localised in passing the metal-insulator transition [139, 140]. 
The relevance of small polaron transport above T c is obvious from the acti- 
vated behaviour of the conductivity [93] . Consequently many theoretical stud- 
ies focused on polaronic approaches [141, 142, 143, 144, 145, 146, 147, 148]. 
Polaronic features have been established by a variety of experiments. For 
example, high-temperature thermopower [149, 150] and Hall mobility mea- 
surements [25] confirmed the polaronic nature of charge carriers in the para- 
magnetic phase. More directly the existence of polarons has been demon- 
strated by atomic pair distribution [151], x-ray and neutron scattering stud- 
ies [152, 153, 154]. Interestingly it seems that the charge carriers partly re- 



Numerical solution of the Holstein polaron problem 



51 



tain their polaronic character well below T c , as proved, e.g., by neutron pair- 
distribution-function analysis [155] and resistivity measurements [156]. 

Regardless of whether the EP coupling acts as a secondary pairing inter- 
action for HTSC in the cuprates, is responsible for the charge ordering in 
the nickelates, or triggers the CMR phenomenon in the manganites, EP and 
particularly polaronic effects need to be reconsidered for the case of strong 
electronic correlations realised in these materials. For instance, Coulomb or 
spin exchange interactions may lead to a "prelocalisation" of the charge carri- 
ers. Then a rather weak EP coupling can cause polaronic band narrowing and 
that way might drive the system further into the strongly correlated regime. 
In the remaining part of this section we will keep track of this problem and 
present some exact results for composite spin/orbital-lattice polarons. 

9.1 Hole polarons in the Holstein t — J model 

Electronic motion in weakly doped Mott insulators like the HTSC cuprates 
is determined by the constraint of no double occupancy of sites and anti- 
ferromagnetic exchange between nearest-neighbour spins. The generic model 
studied in this context is the 2D t — J Hamiltonian, 

H U = -* E + H - c -) + J J2 { S i S j . ( 44 ) 

acting in a projected Hilbert space, i.e. c^J = c-^(l - n it _ a ), n» = J2<? Acfiic^ 
and Si = Y]„„, c'-t.c-,. Within the t — J model the bare transfer am- 
plitude of electrons (t) sets the energy scale for incoherent transport, while 
the Heisenberg interaction (J) allows for spin flips leading to coherent hole 
motion at the bottom of a band with an effective bandwidth determined by 
J. J < t corresponds to the situation in the cuprates, e.g. J/t ~ 0.4 with 
t ~ 0.3 cV is commonly used to model the quasi-2D La2- x Sr x Cu04 system. 

In order to study polaronic effects in systems exhibiting besides strong 
antiferromagnetic exchange a substantial EP coupling the Hamiltionian (44) 
is often supplemented by a Holstein-type interaction term 

H = H tJ ( b l + h i) ~ h i + w ° E ( b l b i + \) (45) 

i i 

(hi = 1 — hi denotes the local density operator of the spinless hole). The 
resulting Holstein t — J model (HtJM) (45) takes the coupling to the hole 
as dominant source of the particle-lattice interaction. In the cuprate context 
an unoccupied site, i.e. a hole, corresponds to a Zhang-Rice singlet (formed 
by Cu3d x 2_ y 2 and 2p XjY hole orbitals) for which the coupling should be 
much stronger than for the occupied (Cu 2+ ) site [157, 158, 159]. The hole- 
phonon coupling constant is denoted by e p = g 2 u>o, and loq is the bare phonon 
frequency of an internal vibrational degree of freedom of lattice site i. 
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The changes of the quasiparticle properties due to the combined effects of 
hole-phonon/magnon correlations are expected to be very complex and as yet 
there exist no well-controlled analytical techniques to address this problem. 
Naturally such a dressed hole quasiparticle will show the characteristics of 
both "lattice" and "magnetic" (spin) polarons. 

We solved the Holstein t — J Hamiltonian on finite lattices using ED, KPM 
and a phonon Hilbert space truncation method. To control our truncation 
procedure we carefully checked the ground-state energy and the weight of the 
m- phonon states (|c m | 2 ) in the ground state as a function of the number of 
phonons retained (M). Convergence is assumed to be achieved if the relative 
error of both E (M) and \c m \ 2 is less than HT 5 . 
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Fig. 29. Phonon-weight function |c m | 2 and ground-state energy Eq(M) for the 2D 
HtJM with J = 0.4 (throughout this section all energies will be measured in units 
of t). 



Figure 29 shows |c m | 2 for the single- hole case at weak, intermediate and 
strong EP couplings. The curves |c m | 2 are bell-shaped and their maxima cor- 
respond to the most probable number of phonon quanta. The importance 
of multi-phonon states becomes apparent especially in the adiabatic strong- 
coupling regime e p 3> t, u>q. 

Let us start the analysis of the 2D HtJM with a discussion of the single- 
hole spectral function 
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Figure 30 displays Ak(u>) for the allowed (nonequivalcnt) momenta K of a 
ten-site square lattice. To visualise the intensities (spectral weights) connected 
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with the various peaks in each iiT-sector we have also shown the integrated 
density of states 



In the absence of EP coupling, of course, we reproduced the single-particle 
spectrum of the pure t — J model [160]. Here one observes a quasiparticlc 
pole corresponding directly to the coherent single-hole ground state (having 
momentum (37r/5,7r/5) on a ten-site lattice) separated by a pseudogap of 
about J from the lower edge of a broad incoherent continuum being ~ 6t 
wide. In the weak coupling regime the mass renormalisation of the coherent 
quasiparticle band due to the hole-phonon Holstein coupling is small compared 
with that arising from hole-spin interactions (magnetic polaron regime). In 
particular, the integrated density of states is barely changed from that of 
the pure t — J model. The new structures, nevertheless observed in the Ak 
spectra shown Fig. 30 (a), correspond to predominantly "phononic" side bands 
separated from the particle-spin excitations by multiples of the bare phonon 
frequency ojq. These phonon resonances have less and less "electronic" spectral 
weight the more phonons are involved. This is because Ak{&) measures the 
overlap of these excited states with the state obtained by creating a hole in 
the zero-phonon Heisenberg ground state. 

With increasing e p the lowest peaks in each if-sector start to separate 
from the rest of the spectrum. These states become very close in energy and 
finally a narrow well-separated lattice hole-polaron band evolves in the strong- 
coupling case (see Fig. 30 (b)). Now the hole is heavily dressed by phonons 
and the quasiparticle pole strength is strongly suppressed (cf. Fig. 31). At the 
same time spectral weight is transferred to the high-energy part and the whole 
spectrum becomes incoherently broadened. Therefore we observe an overall 
smoothing of N(uj). The gap to the next higher energy band is of the order of 
ljo - This excitation will be triggered by an one-phonon absorption process. The 
crossover to the lattice hole-polaron state is accompanied by a strong increase 
in the on-site hole-phonon correlations [38] , indicating that the lattice polaron 
quasiparticlc comprising a self-trapped hole and the phonon cloud is mainly 
confined to a single lattice site (small hole polaron). Most notably, compared to 
the non-interacting single-electron Holstein model [74] or the spinlcss fermion 
Holstcin-t model [161], the critical EP coupling strength for lattice polaron 
formation is considerably reduced due to magnetic prelocalisation effects. 

To illustrate the formation of the lattice hole-polaron band in some more 
detail, we first determined the "coherent" band dispersion Ek of the 2D 
HtJM for a 16-site lattice. The band structure is shown in the left panel of 
Fig. 31 along the principal directions in the Brillouin zone. The minima of 
the quasiparticle dispersion are found to be located at the momenta K = 
(±7r/2, ±7r/2) (the hidden symmetry of the 4x4 cluster leads to an accidental 
degeneracy with the K = (±7r, 0), (0, ±7r) states). At weak EP couplings the 
energy dispersion is not significantly changed from that of the standard t — J 
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Fig. 30. Wavevector resolved single-hole spectral functions Ak{oj) and integrated 
spectral weight N(u>) for the 2D Holstein t — J model at weak (a) and strong (b) 
EP coupling. The results, taken from Ref. [39], were obtained for a tilted \/T0 x 
\/l0 cluster with periodic boundary conditions (if- vectors are given in units of 
(7r/5,jr/5)). 
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model, provided that the phonon frequency exceeds the effective bandwidth 
of the magnetic lattice polaron (lu = 0.8 > AE t j). 

Next we evaluated the wave-function renormalisation factor for different 
band states, 



Z K = 



(^Q,Q^\c i Q-K,<jCQ-K,a\' l f J 0,Q 



(48) 



where \ipQQ ^) denotes the one-hole state being lowest in energy in the Q- 
sector. Zk can be taken as a measure of the "contribution" of the hole (dressed 
at e p — by spin- wave excitations only) to the composite spin/lattice polaron 
(having total momentum K) . The data obtained at weak EP coupling unam- 
biguously confirm the different nature of band states in this regime (see Fig. 31 
right panel). We found practically zero-phonon "hole" states at the band min- 
ima (K = (37r/5, 7r/5), triangles down) and "phonon" states, which are only 
weakly affected by the hole, around the (flat) band maxima (K = (tt,tt), 
triangles up). With increasing e p , a strong "mixing" of holes and phonons 
takes place, whereby both quantum objects completely lose their own iden- 
tity Concomitantly Zk decreases for the "hole-like" states but increases (first 
of all) for the "phonon-like" states. At large e p , a small lattice hole polaron 
is formed, which, according to the numerical data, has an extremely small 
spectral weight. Then the question arises whether the lattice hole polaron is 
a "good" quasiparticle in the sense that one can construct a quasiparticle 
operator, cku — * d,Ka, having large spectral weight at the lowest pole in the 
spectrum. It was demonstrated that such a composite electron/hole-phonon 
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Fig. 31. Left panel: Band dispersion for the 2D HtJM on a 16-site lattice (J = 0.4). 
Right panel: Spectral weight factor Zk as a function of the EP coupling strength 
e p for the different K vectors of the ten-site lattice: (3,1) (triangles down), (2,4) 
(squares), (0,0) (circles), and (5,5) (triangles up) [in units of (7r/5, 7r/5)]. The inset 
shows the narrowing of the coherent bandwidth AE with increasing e p . 
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(polaron) operator could be constructed for the Holstein model [74] as well as 
for the t — J model coupled to buckling/breathing modes [162]. 

In Fig. 32 we show the optical response in the framework of the single- 
hole 2D HtJM. In the weak EP coupling regime and for phonon frequencies 
loq > AE t j, we recover the main features of the optical absorption spectrum 
of the 2D t — J model [163], i.e., an "anomalous" broad mid-infrared band 
(J < u < 2t), separated from the Drude peak (V5{to); not shown) by a 
pseudo-gap ~ J, and an "incoherent" tail up to ui ~ 7t (cf. Fig. 32 (a)). At 
larger EP couplings, the overlap with excited multi-phonon states is enlarged 
and the optical response is enhanced at higher energies. This redistribution 
of spectral weight from low to high energies can be seen in Fig. 32 (b). As 
expected the transition to the lattice hole-polaron state, at about e c v {J = 
0.4, ujo — 0.8) ~ 2.0, is accompanied by the development of a broad maximum 
in a re9 (ui), whereas the Drude weight as well as the low- frequency optical 
response become strongly suppressed. 




Fig. 32. Optical conductivity in the 2D HtJM with J = 0.4. a reg (u) and S re9 (uj) 
were calculated for a ten-site lattice with 15 phonons. 

Let us now make contact with the experimentally observed characteristics 
of the mid-infrared (MIR) spectra in the doped perovskites [22, 160, 164, 165]). 
The simple 2D HtJM seems to contain the key ingredients to describe, at least 
qualitatively, the principal features of the optical absorption spectra of these 
compounds. This can be seen by comparing Figs. 32 (a) and (b) , correspond- 
ing to the weak and strong EP coupling situations realized in the cuprate 
(La2- x Sr x Cu04) and nickelate (La2- x Sr x Ni04) systems, respectively. The EP 
interaction ratio e p /t in the nickelates is estimated to be about one order of 
magnitude larger than in the cuprates because of the much smaller trans- 
fer amplitude (t ~ 0.08 eV [135]). According to the internal structure of the 
low-spin state, the hopping transport of spin-1/2 composite holes in a spin-1 
background is rather complex; implying, within an effective single-band de- 
scription, a strong reduction of the transition matrix elements [166]. A striking 
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feature of the absorption spectra in the cuprate superconductors is the pres- 
ence of a MIR band, centred at about 0.5 eV in lightly doped La2- x Sr x Cu04 
(which, using t ~ 0.3 eV, means that u) ~ 1.5). Such a strong MIR absorption 
is clearly visible in Fig. 32 (a). Obviously it is quite difficult to distinguish 
the spectral weight, produced by the dressing of the hole due to the "bag" of 
reduced antiferromagnetism in its neighbourhood [167], from other (e.g. hole- 
phonon coupling) processes that may contribute to the MIR band observed 
experimentally. The results presented for the HtJM in Fig. 32 (b) support the 
claims, however, that the MIR band in the cuprates has a mainly "electronic" 
origin, i.e., the lattice polaron effects are rather weak [160]. The opposite is 
true for their isostructural counterpart, the nickelate system, where the MIR 
absorption band has been ascribed by many investigators to "polaronic" ori- 
gin [135, 142]. Within the HtJM such a situation can be modelled by the 
parameter set used in Fig. 32 (b). If we fix the energy scale by t = 0.08 eV, 
the maximum in the optical absorption is again located at about 0.5 eV. The 
whole spectrum clearly shows lattice polaron characteristics, where it seems 
that the lattice hole polarons are of small-to- intermediate size [16]. Most no- 
tably, we are able to reproduce the experimentally observed asymmetry in 
the shape of the spectrum, in particular the very gradual decay of a re9 (uj) at 
high energies. It is worth mentioning that this behaviour cannot be obtained 
from a simple fit to the analytical expressions derived for the small polaron 
hopping conductivity [135, 168, 169]. Exploiting the f-sum rule we found that 
there are almost no contributions from band-like carriers in agreement with 
the experimental findings [135, 142]. 

Next let us briefly discuss the two-hole problem. In order to study hole- 
binding effects, we have calculated the hole-hole correlation function 

C ho - ho {\i - j\) = (Me P , J)Millto(e P , J)) • (49) 

Results for C ho _ ho {\i — j\) are presented in Fig. 33. At weak EP coupling, 
Cho-ho{\i — j\) becomes maximum at the largest distance of the ten-site lat- 
tice, while in the intermediate EP coupling regime the preference is on next 
NN pairs. As expected, further increasing e p , the maximum in Cho-ho{\i — j\) 
is shifted to the shortest possible distance, indicating hole-hole attraction. 
At s p » 1, the two holes become "self-trapped" sharing a sizeable com- 
mon lattice distortion, i.e., a nearly immobile hole-bipolaron is formed. The 
behaviour of Cho-ho is found to be qualitatively similar for higher (lower) 
phonon frequencies (see inset), except that the crossings of different hole-hole 
correlation functions occur at larger (smaller) values of e p , which again shows 
the importance of both parameter ratios A = e p /2Dt and g = yJe p /oJo. 

Finally let us consider the quarter-filled band case. Here, we have inves- 
tigated the simpler spinless fermion model (total S z — S z max ). In accordance 
with previous approximate treatments based on an inhomogeneous variational 
Lang-Firsov approach [161], we found, as the EP coupling increases, evidence 
for a transition from a free polaron state to a 2D polaronic superlattice, 
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Fig. 33. Non-equivalent hole-hole pair correlation functions Cho-hoi\i — j\) in the 
two-hole ground state of the HtJM at various e p . Here 1, 2, and 3 label NN, next 
NN, and third NN distances, respectively. 



where the holes are self-trapped on every other site. This crossover is sig- 
nalled by a pronounced peak in the charge structure factor S c (tt, 7t). To visu- 
alise the correlations in this state in more detail, in Fig. 34 we have depicted 
Cho-ho{\i — j\) and the corresponding hole-phonon density correlation func- 
tion Cho-ph(\i — j\) = (V'ol^i^J^'IV'o) as a function of \i — j\. Our exact results 
clearly show the phonon-dressing of the holes and the tendency towards CDW 
formation as observed, e.g., in Lai.sSro.sNiO^ 




Fig. 34. Cho-ho(\6 — j\) (left) and Cho- P h(\6 — j\) (right) are displayed at e p = 3 
and wo = 0.8, where both diameter and gray level of the circles are proportional to 
the correlation strength. 
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9.2 Lattice polarons in a generalised double-exchange model 

The key elements of the electronic structure of the CMR (Lai_ x [Sr, Ca] x Mn03) 
compounds are the partially filled 3d states. The cubic environment of the Mn 
sites within the perovskite lattice results in a crystal field splitting of Mn d- 
orbitals into e g and t 2g . In the case of zero doping (x = 0) there are four 
electrons per Mn site which fill up the three t 2g levels and one e g level, and by 
Hund's rule coupling, form a S = 2 spin state. Doping will remove the electron 
from the e g level, and by hopping via bridging oxygen sites the resulting holes 
acquire mobility. 

Due to the specific symmetry of the manganese d and oxygen p orbitals, the 
transfer of the e g -electrons shows a pronounced orbital anisotropy. In the limit 
of large on-site Coulomb (U) and Hund (Jh) interactions the electron transfer 
is strongly affected by the spin of the core electrons as well. Concentrating 
on the link between magnetic correlations and transport, early studies on 
lanthanum manganites attributed the low-T metallic behaviour to Zener's 
double-exchange (DE) mechanism [170, 171], which maximises the hopping of 
a strongly Hund's rule coupled electron in a polarised spin background. 

As discussed above, it has been realized that physics beyond DE is impor- 
tant not only to explain the phase diagram of the manganites but also the 
CMR transition itself. The orbital degeneracy in the ground state of Mn 3+ 
ions connects the system to the lattice, making it susceptible to Jahn- Teller 
polaron formation. There are two types of lattice distortions which are im- 
portant in manganites. First the partially filled e g states of the Mn 3+ ion are 
Jahn- Teller active, i.e., the system can gain energy from a quadrupolar sym- 
metric elongation of the oxygen octahedra which lifts the e g degeneracy. A 
second possible deformation is an isotropic shrinking of a MnC-6 octahedron. 
This "breathing" -type distortion couples to changes in the e g charge density, 
i.e., is always associated with the presence of an Mn 4+ ion. 

Restricting the electronic Hilbert space to the large Hund's rule states 
given by the spin-2 orbital doublet state 5 E [t 2 ( A A 2 )e] for Mn 3+ (d 4 ) and the 
spin- 1 orbital singlet state 4 A 2 [t 2 ] for Mn 4+ (d 3 ), within 2 nd order pertur- 
bation theory the following Hamiltonian results (for details of the derivation 
and notation see Ref. [172, 173]): 

H = i^DE + ^spin-orbital + -^electron- JT + ^electron- breathing + -ffphonon 




+ ]T (J d KX SiSi+s + AiJP?!*, 
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(50) 



i 



The effective low-energy Hamiltonian H contains Schwinger bosons a\ ' 

i.e. 2Si = a\ fl T llv a i v (fj,,v G {?,!}), fcrmionic holes c^, phonons bfl (a e 

{6, e}), and orbital projectors p^ x *> (k, A e {£, f?, C})- I n Eq. (50), the first 
term, being proportional to t, corresponds to the DE interaction [170, 171]. 
The second term appears to be a bit more involved, since a rather large 
number of accessible virtual excitations (proportional to t 2 and t 2 ) contribute 
(see Refs. [172, 173]). However, in all cases it is basically the product of a 
Hcisenberg-type spin interaction and two orbital projectors. 

The coupling between the orbital degree of freedom of the e g electrons 
and the optical phonon modes to lowest order can be modelled by the E ® e 
Jahn- Teller Hamiltonian (third term) and a Holstein-type interaction (fourth 
term). The energy of the dispersionless optical phonons are given within the 
harmonic approximation (fifth term). 

For analytical methods the above Hamiltonian (50) is far too complex to 
be understood in full detail, and even its numerical solution on finite lattices 
is hard. Using high performance computers and the phonon basis optimisation 
outlined in Sec. 2.2, we were able to calculate the ground-state properties of 
a small four-site cluster and to address, in particular, short-range correlations 
between the charge, spin, orbital and lattice degrees of freedom [172, 173]. 
(See also [62].) Figures 35 and 36 give a glimpse of these results. We assumed 
t = 0.4 cV and t/t n = 3 for the hopping integrals and characterised the 
magnetic "order" according to the total spin of the ground state. 

Undopcd manganites (LaMnC>3, PrMnOs) usually exhibit A- type antifer- 
romagnetic order and strong Jahn- Teller distortion of the ideal perovskitc 
structure. The origin of the observed magnetic order has been subject to dis- 
cussions. While different band structure calculations [174] emphasise the im- 
portance of lattice distortions for the stability of antiferromagnetism, purely 
electronic mechanisms were also favoured [175]. In our microscopic model (50), 
at x = 0, only the second and third term will be active, and without EP in- 
teraction the competition of the spin-orbital contributions depends sensitively 
on the values of Coulomb and Hund's rule coupling. Starting from the "fer- 
romagnetic" phase increasing either U or g changes the magnetic order of 
the ground state to "antiferromagnetism" [172] (see also Fig. 35 (upper left 
panel)). At larger EP interaction the system tends to develop static Jahn- 
Teller distortions, which also fixes the orbital pattern and subsequently the 
spin order. The change of orbital, spin and phonon correlations is illustrated 
schematically in Fig. 36 (left panels). 

Our numerical calculations corroborate the enhancement of ferromagnetic 
correlations for the weakly doped case (see Fig. 35, upper left panel; x — 0.25). 
However, if strong electron-phonon coupling causes self-trapping of the car- 
riers the spin order switches back to antiferromagnetism. This coincidence 



Numerical solution of the Holstein polaron problem 61 

can be seen by comparing the total spin of the cluster and the kinetic en- 
ergy in the ground state, both depicted in the upper panels of Fig. 35. Ob- 
viously the change in the magnetic order is accompanied by the appearance 
of a lattice distortion and a signature in the fluctuation of the bond length 
(oc {Q^iy) — (Qx/y) 2 ), which reminds of the data measured close to the critical 
temperature by Booth et al. [176] (lower left panel). The orbital orientation 
at the sites surrounding the hole-site is sketched in Fig. 36 (middle panels). 
Increasing g isolates the lattice sites, each optimising EP interaction individ- 
ually and uncorrelated with the neighbours. 

U = 6 eV, J = 0.7 eV, t = 0.4 eV, t/t = 3, to. = 70 meV 

n n 
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Fig. 35. Upper panels: Total spin S tot and kinetic energy E kin as a function of 
electron-phonon coupling strength g at various doping levels x. Lower panels: Ex- 
pectation values (q y ) and (q'y) — (q y ) 2 of the bond length in y direction at x = 1/4 
(left) and density-density correlations at x = 1/2 (right). Results obtained by ED for 
the microscopic model (50) on a four site plaquette, where g — g/u>o, and lo — Co = uj 
is assumed [173]. 

At half-filling (x = 0.5) the picture is more involved. Here strong Coulomb 
and EP interactions tend to order the charges in the diagonal direction, i.e., 
in an AB-type structure (compare Fig. 35, lower right panel). This allows 
for a rather large antiferromagnetic spin exchange oc t 2 /Jn- Consequently 
ferromagnetic order is unstable at much lower values of g. The ferromagnetic 
to antiferromagnetic transition is not connected to charge localisation and 
causes only a tiny jump of the kinetic energy. Considering the most relevant 
eigenstate of the bond orbital density matrix, we observe a symmetric order of 
complex orbitals along the diagonal [172]. After charge localisation is achieved 
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Fig. 36. Schematic evolution of charge, spin, orbital and lattice correlation with 
increasing doping [x = 0, x = 0.25 and x = 0.5 (from left to right)] at rather weak 
(upper panels) and strong (lower panels) EP couplings. 

at large g, neighbouring sites are again uncorrelated with respect to orbital 
ordering and are in some real mixed-orbital state. 

Of course, the results of this section will not provide a quantitative anal- 
ysis or description of the real 3D HTSC and CMR materials. Nevertheless, 
ED of even such small systems gives valuable insights into the correlations 
and driving interactions behind the rich phase diagram of the cuprates and 
manganites. Moreover the above ED data may serve as a benchmark for ap- 
proximate theories. 

10 Conclusions 

In summary, we have performed an extensive numerical analysis of the Hol- 
stein model. Combining variational Lanczos diagonalisation, density matrix 
renormalisation group, kernel polynomial expansion, and cluster perturbation 
theory techniques we solved for properties of the Holstein polaron and bipo- 
laron problems. Numerical solution of the Holstein model means that we deter- 
mined the ground-state and low-lying excited states with arbitrary precision in 
the thermodynamic limit for any dimension. Moreover, we calculated the spec- 
tral properties (e.g. photoemission and phonon spectra), optical response and 
thermal transport, as well as the dynamics of polaron formation. Our approach 
takes into account the full quantum dynamics of the electrons and phonons 
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and yields unbiased results for all clcctron-phonon interaction strengths and 
phonon frequencies, but is of particular value in the intermediate-coupling 
regime, where perturbation theories and other analytical techniques fail. 

Most importantly polaron formation represents a continuous crossover of 
the ground state. Nevertheless we found indications for a true phase transi- 
tion in the first excited state, where a polaron plus phonon system changes 
from unbound at weak electron-phonon coupling to bound at strong cou- 
pling. Obviously electron and phonon excitations become intimately, dynami- 
cally connected in the process of polaron formation. Concerning ground-state 
and spectral properties the (quasi-free) electron (small) polaron self-trapping 
crossover is related to (i) a significant increase of the particle effective mass, (ii) 
a substantial narrowing of the polaronic band dispersion, and (iii) a strongly 
suppressed (electronic) quasiparticle residue. At the same time we observe (iv) 
an enhancement of the (on-site) electron-phonon correlations and the forma- 
tion of a phonon drag, (v) a loss of kinetic energy, and (vi) a drop of the Drudc 
weight, accompanied (vii) by a maximum in the spectral weight contained in 
the regular part of the optical response. All these features are found to be 
much more pronounced in higher dimensions. The study of bipolarons showed 
that the two-site bipolaron has a significantly reduced mass and isotope effect 
compared to the on-site bipolaron, and is bound in the strong-coupling regime 
up to twice the on-site Coulomb repulsion. 

Although we have now achieved a rather complete picture of the single 
Holstein polaron and bipolaron problem (perhaps dispersive phonons, longer- 
ranged electron-phonon interaction, finite-temperature and disorder effects de- 
serve closer attention), the situation is disconcerting in the case of a finite car- 
rier density. Here a density-driven crossover from large polarons to quasi-free 
electrons scattered by unbound phonons might occur in the non-adiabatic in- 
termediate coupling regime. Besides electron-phonon coupling competes with 
sometimes strong electronic correlations, e.g., in quasi-lD MX chains, quasi- 
2D high-T c superconductors, 3D charge-ordered nickelates, or bulk colossal 
magneto-resistance manganites. The corresponding microscopic models con- 
tain (extended) Hubbard, Heisenberg or double-exchange terms, and maybe 
also a coupling to orbital degrees of freedom, so that they can hardly be solved 
even numerically with the same precision as the Holstein model. First exact 
results obtained for the case of a single carrier on finite lattices give strong 
evidence that the tendency towards lattice (electron-, hole- or Jahn- Teller-) 
polaron formation is enhanced in strongly correlated electron systems. A more 
thorough investigation of these systems materials and models will definitely 
be a great challenge for solid-state theory in the near future. 
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